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We investigate the dynamics and gravitational- wave (GW) emission in the binary merger of equal- 
mass black holes as obtained from numerical relativity simulations. The simulations were performed 
with an evolution code based on generalized harmonic coordinates developed by Pretorius, and used 
quasi-equilibrium initial data sets constructed by Cook and Pfeiffer. Results from the evolution 
of three sets of initial data are explored in detail, corresponding to different initial separations of 
the black holes, and exhibit between 2-8 GW cycles before coalescence. We find that to a good 
approximation the inspiral phase of the evolution is quasi-circular, followed by a "blurred, quasi- 
circular plunge" lasting for about 1-1.5 GW cycles. After this plunge the GW frequency decouples 
from the orbital frequency, and we define this time to be the start of the merger phase. Roughly 
10-15m separates the time between the beginning of the merger phase and when we are able to 
extract quasi-normal ring-down modes from gravitational waves emitted by the newly formed black 
hole. This suggests that the merger lasts for a correspondingly short amount of time, approximately 
0.5-0.75 of a full GW cycle. We present first-order comparisons between analytical models of the 
various stages of the merger and the numerical results — more detailed and accurate comparisons will 
need to await numerical simulations with higher accuracy, better control of systemic errors (including 
coordinate artifacts), and initial configurations where the binaries are further separated. During the 
inspiral, we find that if the orbital phase is well modeled, the leading order Newtonian quadrupole 
formula is able to match both the amplitude and phase of the numerical GW quite accurately until 
close to the point of merger. We provide comparisons between the numerical results and analytical 
predictions based on the adiabatic post-Newtonian (PN) and non-adiabatic resummed-PN models 
(effective-one-body and Pade models). For all models considered, 3PN and 3.5PN orders match the 
mspiral numerical data the best. From the ring-down portion of the GW we extract the fundamental 
quasi-normal mode and several of the overtones. Finally, we estimate the optimal signal-to-noise 
ratio for typical binaries detectable by GW experiments. We find that when the merger and ring- 
down phases are included, binaries with total mass larger than AQMq (sources for ground-based 
detectors) are brought in band and can be detected with signal-to-noise up to ~ 15 at 100 Mpc, 
whereas for binaries with total mass larger than 2 x l(f Mq (sources for space-based detectors) the 
SNR can be ^ lO'' at IGpc. 

PACS numbers: 04.25.Dm, 04.30.Db, 04.70.Bw, 04.25.Nx, 04.30.-w 

I. INTRODUCTION 

With gravitational- wave (GW) detectors operating or under commissioning |4j], it is more and more desirable 

to improve the theoretical predictions of the GW signals. Compact binaries composed of black holes (BH) and/or 
neutron stars are among the most promising candidates for the first detection. 

The last year has been marked by breakthroughs in numerical relativity (NR) with several independent groups 
being able to simulate binary black hole coalescence through the last stages of inspiral (2-4 orbits), merger, ring-down 
and sufficiently long afterwards to extract the emitted GW signal 0, 0, 0, d, Q . In Ref . Q the first stable evolution of 
such an entire merger process was presented. A couple of the key elements responsible for this success where the use 
of a formulation of the field equations based on a generalization of harmonic coordinates [13, [HI and the addition 
of constraint damping terms to the equations ^13. , 14]. Similar techniques have been successfully incorporated into 
other efforts since [TEl [TgI [TtI . Some months afterwords, two groups [a 0] independently presented modifications of 
the BSSN (or NOK) [l8l. Il9l. [20l| formulation of the field equations that allowed them to simulate complete merger 
events. Among the key modifications were gauge conditions allowing the BHs to move through the computational 
domain in a so-called puncture evolution [21]; these methods have also been successfully reproduced by other groups 
since BilHiilll. 

In this paper, after analyzing the main features of the dynamics and waveforms obtained from the NR simulations, 
we present a preliminary comparison between the numerical and analytical waveforms. The analytical model we use 
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for the inspiral phase is the post-Newtonian (PN) approximation [25|, |26[ in the adiabatic hmit ^ [26[ , whereas for 
the inspiral-(plung;e)-mer ger-ring -do wn p hases we consider the PN non- adiabatic and resummed models, such as the 
effective-one-body (EOB) |27l. l28l [29l. ISOl. Isil Is^] and Pade resummations [13] ■ Due to the Umited resolution, initial 
eccentricity, relatively close initial configurations and possible coordinate artifacts, it is difficult to claim very high 
accuracy when comparing with analytical models. Thus, we shall refer to those comparisons as first-order comparisons. 
We will only consider a few dynamical quantities of the analytical models characterizing the binary evolution, notably 
gauge-invariant quantities such as the orbital and wave frequencies, and the GW phase, all as measured by an observer 
at infinity. When comparing, we will assume that the numerical and analytical waveforms refer to equal-mass binaries, 
but we apply a fitting procedure to obtain the best-match time of coalescence and the spin variables. We will use the 
confrontation with analytical models as an interesting diagnostic of the numerical results. When simulations that are 
more accurate and begin closer to an inspiralling circular binary become available we will be able to do more stringent 
tests of analytical models, compare all dynamical quantities expressed in the same gauge, and use those results to 
discriminate between models. 

Certainly, the most intriguing and long-awaited result of the numerical simulation is the transition inspiral-merger- 
ring-down. Is it a strongly non-linear phase? How much energy and angular momentum is released? Over how many 
GW cycles does it occur? How spread in frequency is the signal power spectrum? Answers to these questions are 
relevant from a theoretical point of view, e.g., to study general relativity in the strongly coupled regime, and also 
from an observational point of view, e.g., to build faithful templates to detect GW waves and test GR with GW 
experiments. In this paper we shall start to scratch the surface of this problem. We pinpoint several interesting 
features of the inspiral to ringdown transition that need to be investigated more quantitatively in the future when 
more accurate simulations become available. 

Quasi-normal modes (QNM) (or ring-down modes) of Schwarzschild and Kerr BHs were predicted a long time 
ago [3^, [13, nil . Their associated signal can be described analytically in terms of damped sinusoids. By fitting 
to the numerical waveforms we extract the dominant QNMs of the final Kerr BH — i.e. the fundamental mode and 
several of the overtones — and try to make connections to the previous dynamical phase. 

Finally, we discuss the impact of the merger and ring-down phases on the detectability of GWs emitted by equal- 
mass binaries for ground-based and space-based detectors, and compare those results with predictions from analytical 
models. 

This paper is organized as follows. In Sec. |TT] we review the initial-data sets of Cook and Pfeiffer [Zoj used in the 
numerical simulations. In Sec. lIIII we present and discuss the results of NR simulations of binary BH mergers obtained 
with the generalized-harmonic-gauge code of Pretorius Q. In Sec. IIVI we provide a first-order comparison between 
numerical and analytical results for the last stages of inspiral. In Sees. |V] and IVII we analyze the ring-down and 
merger phases as predicted by the numerical simulations. In Sec. IVIll we present a first-order comparison between the 
numerical results and the EOB predictions for inspiral, plunge, merger and ring-down. In Sec. IVIIII we evaluate the 
Fourier transform of the waveforms and discuss how the inclusion of the merger and ring-down phases will increase 
the optimal signal-to-noise ratio of ground-based and space-based detectors. Section |IX] contains our main conclusions 
and a discussion on how to make more robust comparisons with analytical models in future NR simulations. 

Some material we defer to appendices. The majority of the comparisons and analysis of gravitational waveforms 
focus on the dominant quadrupole multipole moment of the wave; in Appendix|^we briefly describe the sub-dominant 
multipole moments extracted from the waves. Appendix IB] compares the energy and angular momentum flux of the 
numerically extracted GW with analytical models of the fluxes. In Appendix [Cl we describe some possible artifacts 
induced by extracting the GW a finite distance from the source in the simulations. Appendix ID] contains tables of 
fitting coefficients from the QNM ring-down fits. We have analyzed three sets of BBH evolutions — some figures from 
the case with the closest initial separation are contained in Appendix |E] to simplify the main text. 

II. INITIAL DATA 

The evolutions presented in this pap er begin with initial data that has been prepared using the methods developed 
by Cook and Pfeiffer [iO, ST], |4^. l43j. This approach incorporates the extended conformal thin-sandwich (CTS) 
decomposition [S, '45*1 , the "Komar mass method" for locating circular orbits [i^, H^l , and quasiequilibrium boundary 
conditions on BH excision surfaces ^40^ A9\ . The data we have used are designed to represent an equal- mass binary BH 
configuration in which the binary is in quasiequilibrium with the holes in a nearly circular orbit and where the spins 
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of the individual holes correspond to "corotation" . Within the CTS approach, the conformal metric, the trace of the 
extrinsic curvature, and their time derivatives must be freely specified. The time derivatives of the conformal metric 
and the trace of the extrinsic curvature are chosen to vanish. This time derivative is taken along an approximate 
helical Killing vector which defines the notion of quasiequilibrium. The initial data is constructed on a "maximal 
slice" which fixes the trace of the extrinsic curvature to zero. Finally, the conformal metric is chosen to be flat. 

The initial data produced by this procedure do a very good job of representing the desired astrophysical situation 
of a pair of BHs nearing the point of coalescence. However, two important approximations have been made in 
the construction of these initial-data sets. When the initial data are evolved, these approximations will affect the 
subsequent dynamics and the GW that is produced. The first approximation is that the initial data are "conformally 
flat" . The choice of a flat conformal metric is known to introduce small errors in representing both individual spinning 
BHs and binary systems. As part of this error, some amount of unphysical gravitational radiation is included in the 
initial data. The second approximation is in placing the binary in a circular orbit. This approximation is motivated by 
the fact that for large enough separation, the time scale for radial motion due to radiation reaction is large compared 
to the orbital period. However, this approximation results in BHs having little, if any, initial radial momentum. For 
sufficiently small separations, this is clearly not "astrophysically correct." 

Until now, the best way to estimate the quality of BH binary initial data has been to compare them against the results 
obtained by FN methods. Comparisons of gauge- invariant quantities such as the total energy and angular momentum 
of the system or the orbital angular velocity, all measured at infinity, are in good agreement with adiabatic sequences 
of circular orbits as determined by third-order FN (3FN) calculations [l^ and their BOB resummed extension [23, H^l 
(see Figs. 10-18 in Ref. [13], Figs. 3-5 in Ref. [131 and Figs. 3-5 in Ref. [5l| and discussion around them). There 
are, of course, differences between the numerical initial data and the PN/EOB models and these differences increase 
as the binary separation decreases and the system becomes more relativistic. Among these differences, there is also 
some evidence that the initial orbit incorporates a small eccentricity (s^ . 

Adiabatic sequences of BH binaries exhibit an inner-most stable circular orbit (ISCO) defined by a turning point in 
the total conserved energy. For the numerical initial data, the quasiequilibrium approximation becomes less accurate 
as the binary separation decreases and we would expect the approximation to be rather poor at the ISCO. FN methods 
restricted to circular orbits suffer similar problems as they approach the ISCO. The PN/EOB and numerical initial- 
data circular-orbit models are in reasonable agreement up to the ISCO, but it is difficult to ascertain the accuracy of 
either in this limit. 

The comparisons done in Refs. [13, H^i HII between numerical initial data and PN/EOB models are limited by 
the fact that they use adiabatic circular orbits. Essentially, these comparisons lend strength to the belief that the 
conformal-flatness approximation is not causing significant problems with the non-radiative aspects of the initial 
data. Clearly, they cannot shed any light on the effects of the circular-orbit approximation. Some insight on this 
issue has been obtained by performing full dynamical evolutions of the PN/EOB equations of motion for equal-mass 
binaries (28l. [ssj. These studies show that neglecting the radial momentum, at both the initial time and throughout 
the evolution as done in adiabatic circular orbits, results in a phase error in the waveforms (see e.g.. Fig. 5 in Ref. [28j). 
Neglecting the radial momentum at the initial time also introduces eccentricity into the dynamics (see e.g., Fig. 4 in 
Ref. [13). 



III. NUMERICAL RELATIVITY RESULTS AND THEIR DIAGNOSTICS 



A. Initial Data for Generalized Harmonic Evolution 



The corotating quasi-circular BH inspiral data discussed in the previous section were evolved using a numerical code 
based on a generalized harmonic (GH) decomposition of the field equations, as described in detail in Ref. [sl [s^. [ssj. 
As supplied by Pfeiffer 42, 43], the initial data is given in terms of standard 3 -I- 1 or ADM (Arnowitt-Deser-Misner) 
variables [56l l57l. ISSi. [S^. namely the lapse function a, shift vector spatial metric hij and extrinsic curvature Kij: 

ds^ = -a^dt^ + h,j {dx' + P'dt) {dx' + 13' dt) , (1) 
= -h/h,"'V^ni. (2) 

In the above — ~aSI f^t is the unit time-like vector normal to t — const, hypersurfaces, and we use units where 
G = c = 1. The GH code directly integrates the 4-metric elements g^i, 

ds^ = g ^^dx^ dx'^ , (3) 

and therefore needs the values of g^^ and dgfi^/dt at i = as initial conditions. The initial data ([T]), ^ provides 
most of what is required to construct g^v\t=OTdg^v/dt\t=Q] what must still be specified are the components of the 
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gauge encoded in da/dt\t=[j and 8(3^ /dt\t=o. We choose the time derivatives of the lapse and shift such that the shce 
is spacetime harmonic at t = 0: 

d^a 11? = —aK (4) 

d^P' = afikh''' - djU (5) 

where K is the trace of the extrinsic curvature, and F^ j. are the Christoffel symbols of the spatial metric hij. 



B. Characterization of the Waveform 



Gravitational wave information is obtained by computing the Weyl scalar ^'4, which has the asymptotic property 
of being equal to the outgoing radiation if the complex null tetrad is chosen correctly. To be explicit, we define 
a spherical coordinate system centered on the center of mass of the binary with orthonormal bases (r,0, 0). The 
coordinates are chosen so that the azimuthal axis is aligned with the orbital angular momentum and the binary orbits 
in the direction of increasing azimuthal coordinate. 

To define our complex null tetrad, we use the time-like unit vector normal to a given hypersurface fi and the radial 
unit vector f to define an ingoing (fc) and outgoing null vector {€} by 

k = -^('^ + ^)^ (6) 
i=(n-f). (7) 



We define the complex null vector rfi by 



In terms of this tetrad, we define ^'4 as 



i=(^_^^). (8) 



^-4 = Ca.f3jse"{m^re-'{m'r, (9) 



where Capjs is the Weyl tensor and * denotes complex conjugation. 

To relate ^'4 to the GWs, we note that in transverse-traceless (TT) gauge, 



Following convention, we take the and polarizations of the GW to be given by 

^+ - ^(^ff-^JJ)^ (12) 
h. = hJl (13) 

We find, then, that in vacuum regions of the spacetime, 

VI/4 = /i+ - i/ix. (14) 

It is most convenient to deal with ^4 in terms of its harmonic decomposition. Given the definition of ^'4 in Eq. ^ 
and the fact that to* carries a spin- weight of —1, it is appropriate to decompose \l/4 in terms of spin- weight —2 
spherical harmonics _2y£m(0, 0). There is some freedom in the definition of the spin-weighted spherical harmonics. 
To be explicit, we defined the general spin- weighted spherical harmonics by 



/o/ 1 1 

sYe^{0,4,) EE {-iyJ—^dl^_,^^{e)e"'^f, (15) 
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where d^^^ is the Wigner d-function 



(~i) + m)\{i - m)\{£ + sy.je - sy . ^ 

^ {£ + m- ty.{e -s- t)\t\{t + s-m)\ ' 



t=Ci 



and where Ci — max(0, m ~ s) and C2 = min(£ + m, ^ — s). 

Finally, for convenience, we always decompose the dimensionless Weyl scalar where Al = mi + m2 is the 

mass of the initial binary system with mi and m2 the irreducible |60j] masses of the individual BHs, and r is the 
generalized harmonic radial coordinate. We then define 

rM*4(i,r) = ^ -2C,™(t, r)_2y,™(0, 0). (17) 

The complex mode amplitudes -~2Cernit, r), extracted at a fixed generalized harmonic coordinate radius r, contain the 
full information about the gravitational waveforms as a time-series. 

In the numerical code the four orthonormal vectors {h, f, 9, (f) used to construct the null tetrad are computed as fol- 
lows. The spacetime is evolved using Cartesian coordinates x, y, z with time and we use the standard transformation 
to define the spherical coordinates: 

X — r cos((/)) sin(6') , (18) 
y = r sin(0) sin(0) , (19) 
z = rcos{e), (20) 

n is the time-like unit vector normal to t = const, surfaces, and f is the unit space-like vector pointing in the direction 
[d/drY. In the limit r — > cx) the time coordinate t coincides with TT time. 9 is computed by making (9/(90)" 
orthonormal to (n,f) using a Gramm-Schmidt process, and then (p is calculated by making {d/d(j))°' orthonormal to 
(n, f, 9). All norms are computed with the full spacetime metric 5^1,. The Weyl scalar ^1*4 is evaluated over the entire 
numerical grid (i.e. at all x, y, z mesh points) at regular intervals in time. We then interpolate ^'4(t, x, y, z) to a set of 
coordinate spheres at several "extraction radii" r^, with a uniform distribution of points in (9, (f>) ^. All the waveform 
related data from the simulations presented here are taken from such samplings of ^4(<, r — ri,9^(j)), and we have 
used ri = 12.5m, 25m, 37.5to and 50to. The plots and comparisons shown in the main part of the paper use = 50m, 
while in Appendix [Cl we discuss the trends that are seen in '^4 as a function of r^. To summarize the results of the 
Appendix, prior to merger the extraction radii at 37.5m and 50m appear to be well within the "wave-zone" , and thus 
gives a decent representation of the waveform. Specifically, the coordinate propagation speed of the wave from one 
extraction sphere to the next is very close to unity, and the structure of the wave, normalized by r, is similar at the 
two extraction points. Interestingly though, at later times during the simulation, apparently in coincidence with the 
strongest wave emission around the merger part of the evolution, the gauge seems to change slightly in the extraction 
zone. The coordinate speed drops by several percent, and the amplitude of the normalized waveform also decreases 
as the wave moves outward. The effect is more pronounced for binaries that are initially further separated. As not 
all metric information was saved when the simulations were run we cannot describe what underlying properties of 
the metric are responsible for the change in propagation behavior, though for the purposes of this paper the effect is 
sufficiently small that we do not believe it will alter any of the primary conclusions. A couple of exceptions are in 
the estimates of the total energy and angular momentum radiated — these calculations involve double and triple time 
integrations of squares of the wave, and so significantly amplify even small systematic errors in the waveform. This 
is, we believe, the cause of the increasing over-estimate of these quantities as the initial orbital separation increases, 
as shown later in Table HIl Future work will attempt to address these gauge- related issues. 



Numerical results and a discussion of errors 



We have evolved 3 sets of initial data, labeled hy d = 13,16 and 19 in Ref. [401; the initial orbital parameters 
are summarized in Table U'^. Each initial data set was evolved using three different grid resolutions, summarized 



33 X 65 points in 9 £ 0..7r,</) G 0..27r for this set of simulations 
^ Note that when writing 81^2/^1 2 rnean by mi 2 the irreducible mass nij^'j- When later on we compare with PN models, we should 
first express the rest masses appearing in the PN formula in terms of m'j^j s-n^i the spin variable 5i,2- However, since we deal here with 
very small spins, the error in not doing so is very small. 
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"d" 


Madm/M 


Mlu 


Jadm/M^ 




Si /mi 


13 


0.986 


0.0562 


0.875 


7.96 


0.107 


16 


0.988 


0.0416 


0.911 


9.77 


0.0802 


19 


0.989 


0.0325 


0.951 


11.5 


0.0629 



TABLE L Several parameters describing the initial data [iol ] (left to right): the ADM mass of the spacetime, the initial angular 
velocity of each BH, the ADM angular momentum, the initial proper separation between the holes, and their initial spins. The 
units have been scaled to AI, the sum of initial AH masses. 



"Resolution" 


wave-zone res. 


orbital-zone res. 


BH res. 


h 


1.6Af 


0.20M 


0.048M 


3/4 h 


1.2Af 


0.15M 


0.036M 


1/2 h 


0.82M 


O.lOAf 


0.024M 



TABLE II: The three sets of characteristic spatial resolutions used in the simulation discussed here, where each resolution is 
labeled relative to the coarsest resolution h. The grid is adaptive with a total of 8 levels of refinement, and the coordinate 
system is compactified. The wave zone is defined to be at r = 50M, the orbital zone within about r — lOM and the BH zone 
is within 2 — 3A/ of each AH. A Courant-Friedrichs-Lewy (CFL) factor of 0.2 was used in all cases. 



in Table [TTl Most of the results presented in this paper are from the highest resolution simulations, with the lower 
resolution runs providing error estimates via the Richardson expansion. During evolution the same temporal source- 
function evolution equations were used as with the scalar-field-coUapsc binaries described in Ref. [5^ . and the 
(covariant) spatial source-functions were kept equal to zero. We have quite extensively tested this code to make 
sure we are solving the Einstein equations, and convergence tests of the residual of the Einstein equation for similar 



simulations were presented in Ref. 55 | 



Table Hill lists some key information obtained from each merger simulation. Figure [T] shows the orbital motion prior 
to merger. The uncertainties and error estimates listed in the table were calculated using an assumed Richardson 
expansion, as discussed in '55||. At least three simulations are needed to verify that one is in the convergent regime 
in general, and for the properties listed in the table we do see close to second order convergence for most properties 
— a couple of anomalous cases (the phase error and merger time for the d = 19 case) are discussed further below. 
Assuming one is in the convergent regime, results from two simulations can then be used to estimate the truncation 
error. This estimate will not account for systematic errors, and here we list several potential sources of such error. 
Some quantities are in principle susceptible to gauge or coordinate effects, including the apparent horizons (AHs), 
and therefore properties measured using them; orbital parameters such as angular frequency and eccentricity deduced 
from the positions of the AHs; the finite GW extraction radius and nature of the coordinates at the extraction surface 
(see Sec lIII Bl and AppendijCj): and the choice of tetrad used to calculate \l/4 [ill, IH, HI] . Certain post-processing 
operations have a truncation error associated with them and we have not estimated the magnitude of these. For 
example, in all cases ^4 was sampled on a mesh of size 33 x 65 points in x </> at the extraction radius at a given time. 
We suspect that many of these systematic errors are small, though eventually the accuracy of the simulations will be 
improved (through a combination of high-order methods and higher resolution), and then it may become important 
to quantify and eliminate these additional uncertainties. 

In a waveform, the error of most significance to GW detection is an error in the phase. The first generation of 
numerical binary BH merger simulations [i, H, 0, B m, m, m, m , with the notable exception of the Caltech/Cornell 
effort (66j . all suffer from rather significant cumulative phase errors in the inspiral portion of the waveform for 
the longer duration merger events. In Ref. [11] it was argued that the dominant portion of the phase could be 
factored out as a constant phase shift within the wave, resulting in a "universal" merger result for the class of initial 
conditions considered. Henceforth, we shall denote by maquillage the operation of phase shifting waveforms to make 
them agree at a specified point in time, improving their coincidence (appeareance) over a larger interval of time. 
For certain applications this maquillage waveform is the relevant one. For example, in a matched-filter search the 
initial phase of the waveform is an extrinsic parameter [3^ and is irrelevant for detectability of the signal. Also, 
when comparing waveforms from different initial-data sets the waveforms need to be aligned in some manner for a 
meaningful comparison, and so again a constant phase difference, whatever the source, is largely irrelevant. However, 
for parameter estimation in an inspiral search where a hybrid PN/numerical template is used, the association of 
a numerical merger/ring-down to a FN inspiral waveform will be very sensitive to all phase errors. In particular, 
an uncertainty in the overall phase evolution prior to merger in the numerical waveform is directly related to an 
uncertainty in the merger time for the given initial conditions, and this will translate to an uncertainty in the FN 
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d=13 


d=16 


d=19 


Mf/M 


0.950 ± 0.005 


0.954 ± 0.005 


0.952 ± 0.005 


Egw/M 


0.036 ± 0.004 


0.043 ± 0.004 


0.052 ± 0.004 


af/Mf 


0.71 ±0.02 


0.71 ±0.02 


0.70 ± 0.02 


J}/M^ 


0.64 ±0.02 


0.65 ±0.02 


0.63 ± 0.02 




0.23 ±0.02 


0.31 ±0.02 


0.42 ±0.02 


number of orbits 


1.47 ±0.10 


2.47 ±0.09 


4.39 ±0.18 


tcAu/M 


109 ±4 


228 ± 16 


529 ± 22 


(ipeak — tcAtl)/M 


^ 9 


f» 9 


f» 9 


(idee - tcAH)/M 


^ -11 


^ -11 


fa -9 


initial eccentricity ei 


— 


— 


0.018 ±0.003 


initial eccentricity 62 






0.012(±0.014, -0.012) 


Max. GW amp. error 


8% 


9% 


8% 


(Max. GW phase error)/(27r) 


0.08 


0.7 


1 


(Max. "shifted" GW phase error) /(27r) 


0.04 


0.06 


0.05 



TABLE III: A summary of simulation results. From top to bottom: (i) the mass Alf of the final BH estimated from AH 
properties, (ii) the energy Egw emitted in GWs extracted using the Newman-Penrose scalar ^'4 at a coordinate radius of 
r = 50M from the origin, (iii) the Kerr spin parameter a/ and corresponding non-zero z-component of the angular- 
momentum vector of the final BH (from AH properties), (iv) an estimate Jew of the angular- momentum radiated in GWs 
(from *I'4), (v) the number of orbits in coordinate space before a common AH, at coordinate time t = tcAu, is first detected, 
(vi) an estimate t = tpcak of when the GW amplitude reaches its peak, (vii) an estimate t — tdcc of when the GW frequency 
decouples from the orbital frequency, (viii) for d = 19 two estimates of the eccentricity in the initial data calculated using 
(|21|) (ei) and (|22p (62), (rx) the estimated maximum error in the amplitude of the extracted GWs (occurring near the peak of 
emission), (x) an estimate of the cumulative phase error in the wave from t = until merger, and (xi) the phase error after 
maquillage. The estimated uncertainties and errors do not include possible systematic errors; see the discussion in Sees. IIII Bl 
nil Gl and Appendix [C] for more information. In all plots in the paper where we have shifted time by either tcAH or tpeak, any 
quantity shown that was measured from the waveform at r = 50M is also shifted by the propagation time 5t of the wave to 
r = 50M. We have estimated St by finding the shift resulting in a best-fit between orbital frequency measured using the GW 
versus AH orbital motion, as shown in Fig. [T] specifically, we get 5t = 65M, 68Af and 70Af for d — 13, 16 and 19 respectively. 

binary parameters identified with the match. 

For the reasons just outlined, in Table IIIII we give two estimates of the phase errors in the waveforms. The first 
is the cumulative error in the phase directly measured from the waveform, and the second is the cumulative error 
after the maquillage. Specifically, in the latter case we shift all the waveforms in time so that the peak amplitudes 
(corresponding to the peak of the energy radiated) occur at t=0, and then apply a constant rotation in the complex 
plane of the waveform to give optimal overlap with a reference waveform (typically the highest resolution result). 
Example waveforms before and after the shift for the three resolutions are given in Figs. O [31 Unfortunately, the 
lowest resolution waveform data for the d = 19 case was accidentally deleted, and so only the medium and higher 
resolution results are shown. The two d = 19 cases are unusual in that the phase difference between them does not 
grow monotonically with time; rather the phase difference initially grows then decreases so that by merger time the 
difference is close to zero. We are not entirely sure why the error in the phase evolution behaves so in this case. 
One possibility is that for this longer integration time we are too far from the convergent regime to see the trends in 
phase evolution observed for the shorter d=13 and d=16 runs. However, other estimators of convergence, including 
AH properties and orbital motion as shown in Figs. 31 [S] which do include data from the lowest resolution suggest the 
d = 19 case is in the convergent regime. Another possibility is that the numerical error in the phase has a periodic 
time component whose frequency is sufficiently low that the d=13 and 16 runs do not show it. Regardless, that the 
l/2h and 3/ Ah d=19 cases merge at almost the same time makes it impossible to use them to estimate the error in 
merger time and phase. Therefore, the errors quoted for these two numbers in Table Hill for d=19 were estimated 
using the difference in merger time between the h and l/2h d=19 runs, with the estimated phase error being the error 
in the merger time divided by the waveform period at ring-down multiplied by 27r. 

A final note regarding the phase difference in maquillage waveforms: for two orbits that are quite similar, ei- 
ther because of similar initial conditions or the same initial conditions but different numerical truncation error, the 
time/phase shifting can be performed at any time during the common inspiral/merger phase. The phase difference 
will then by construction be identically zero at the time of the match, and slowly drift as one moves away from the 
matching time. The choice of matching at the peak of the wave amplitude effectively minimizes the net phase error 
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FIG. 1: (left) The orbital motion of one BH from each of the three cases: the dot-dashed (blue) line is d = 13, the solid (black) 
d = 16 and the dashed (red) d — 19. The position of the BH is defined as the center of its AH, and the curve ends once an 
encompassing horizon is found. The eccentricity present in the initial data is particularly evident for the d = 19 case, though 
in part this is due to numerical error — see Fig. [5] To aid in the comparison each trajectory was rotated by a constant phase so 
that they coincide at a coordinate separation of 3M. (right) The orbital motion of the BH's from the d=16 simulation, showing 
the coordinate shapes of the AHs at several key moments. Also shown is the location of the co-rotating light ring of the final 
BH — see the discussion in Sec. IVII 



as this is where the wave frequency is highest. 



D. Diagnostic of the orbital evolution 



The initial orbits of the d — 19 case displayed in Fig. [T] are clearly neither circular nor a smooth adiabatic inspiral. 
It is natural to refer to such orbits as being eccentric. However, describing orbits as "eccentric" when radiative effects 
are strong can be problematic. The notion of eccentricity is precise in Newtonian physics, where the eccentricity is one 
of two parameters needed to describe a general, bound elliptic orbit. In general relativity, even when considering only 
the conservative dynamics, binaries do not follow closed elliptic orbits. When the dissipative effects of gravitational 
radiation are strong, it becomes even more difficult to define the concept of eccentricity. 

The initial data we use starts with essentially no radial momentum. If radiative dissipation is neglected, such orbits 
can be circular or eccentric depending on the magnitude of the orbital angular velocity. But, because of radiation 
reaction, initial data with no radial momenum cannot represent a binary on a smooth quasi-circular inspiral. In fact, 
such an orbit must have some effective eccentricity. 

We will use two methods to attempt to calculate this eccentricity in the d — 19 case. Neither of these methods 
work well for the d = 13 and rf = 16 cases as they do not exhibit enough orbital motion prior to merger. The first 
method uses the following relationship that holds for an orbit with eccentricity e, orbital angular frequency to and 
separation r in Newtonian theory: 

Lj^{t)r^{t)/M - 1 = ecos(0(i)). (21) 

Lower order general relativistic corrections (in particular perihelion precession) will change the argument to the 
cosine function, though the amplitude remains e. The right panel in Fig. [5] shows the LHS of Eq. (|2ip for the d = 19 
simulation, with oj and r calculated from the coordinate motion of the BHs. A certain amount of eccentricity is due 
to numerical error, though the trend in the curves of Fig. [5] as resolution increases indicates that some amount of 
eccentricity does come from the initial data. The numerical data does not follow Eq. (pij) too closely, though at early 
times there are clear oscillations about a line, and we will use the amplitude of these oscillations to define e. For a 
fitting function we use ao + ait + e cos{a2t + 03), and guided by Eq. (|2ip. we define the amplitude of the oscillation to 
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FIG. 2: The Re[_2C22] component of the d = 13 and d — 16 waveforms, unshifted (top) and shifted (bottom). AU resolutions 
are shown to demonstrate the size of numerical errors in the simulation, and data such as this was used to compute the errors 
for waveform quantities listed in Table HIl 
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FIG. 3: In the left panel we show the the Re[_2C22] component of the d = 19 waveform, unshifted (top) and shifted (bottom). 
The lower resolution data for the d = 19 case was accidentally deleted. In the right panel we zoom in for a close-up of the inspiral 
part of the unshifted waveform. From these two results alone it would appear as if the d = 19 simulations have anomalously 
good convergence behavior (compare to Fig. [Sjl. However, this is not the case — refer to the discussion in Sec. IIII Cl and see 
Figs. |4] and [5] for other estimators of the convergence of the solution. 
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FIG. 4: Sum of AH masses (left panel), and the Kerr angular momentum parameter of the final BH (right panel) for the 
d = 19 simulations. The angular momentum was estimated using the ratio of polar to equatorial proper circumference of the 
horizon [63]; the dynamical horizons estimate [o^l gives similar results modulo the oscillations about the mean. Except near the 
time of merger the sum of AH masses in the spacetime should be conserved, and similarly at late times for the Kerr parameter. 
As resolution increases we see the expected trends in these quantities. Note that the "jaggies" in the AH mass estimate is a 
refiection of AH finder problems in the code, and not irregularities in the underlying solution [ssj . 




FIG. 5: In the left panel we show the coordinate separation of the BHs as a function of time for the d = 19 simulations. This 
plot highlights the eccentricity within the orbit and it also refiects the phasing behavior in the waveform for the 3/4/i and l/2h 
cases — see Fig. [S] In the right panel we estimate the eccentricity for the d = 19 case: shown is a plot of the left hand side of 
Eq. (|2H) together with a fit of the form ao + ait + e cos{a2t + as) to the early time behavior of this function. We estimate the 
eccentricity to be the amplitude of the sinusoidal part of the fitting function. 



11 




-500 -400 -300 -200 -100 
(t - V, J/M 



FIG. 6: We plot the ratio between the radial velocity and the tangential velocity as evaluated in the numerical evolution of the 
d = 19 case (dashed curve), and compare it with the Newtonian predictions (dot and continuous curves) obtained using the 
quadrupole formula. The two dashed vertical lines span the region in which a dynamical ISCO could be present. 



be the eccentricity. For the l/2h run the fit gives e = 0.018 ± 0.003, with the uncertainty calculated using the l/2/i 
and 3/4/i data and assumed second order convergence. 

For a second estimate of the eccentricity we use another Newtonian definition given in Ref. ^69] : 

/77i — ATi 

(22) 



where ujp is the frequency w at a local maxima, and LUa is the frequency at the following local minima. Using this 
definition, and the data for the d = 19 case shown in Fig. [3 we get e = 0.012 (0.029,0.068) for the l/2h (3/4/i, h) 
resolution runs. The rather large differences in the values calculated using the different resolutions means that the 
corresponding uncertainty in e calculated using Eq. is also large: ±0.014 (for the values quoted in Table Hill we 
restricted e > 0). 

As mentioned in Sec. |TT1 Refs. [s^ have shown that 3PN estimates of eccentric orbits suggest the quasicircular 
initial data being used has some intrinsic eccentricity. From Fig. 2 of Ref. [s^], we find a 3PN estimate of e ~ 0.01 for 
the d = 19 case. We note that this is remarkably close to the eccentricity estimate obtained via Eq. ((^^ . However, 
despite this coincidence, we should be cautious in attributing the "eccentricity" observed in the orbit of the d = 19 case 
to a non-vanishing eccentricity in the initial data. It is important to remember that the initial data are constructed 
to have vanishing radial velocity. As shown by Miller [s^ , initially circular orbits clearly lead to the kind of effective 
eccentric behavior seen in our numerical evolutions. A comparison of Fig. 4 of Ref. [H^l with Fig. Oof this paper also 
shows striking similarity. It is clear that the initial data, through a combination of a vanishing initial radial velocity 
and possibly non- vanishing initial eccentricity, results in an evolution that exhibits some undesired eccentric behavior. 
However, it is not yet possible to determine which, if either, effect dominates. 

Despite the presence of the eccentricity, the orbital motion on average is quasi-circular. By this we mean that 
throughout the evolution the radial velocity is smaller than the tangential velocity. At leading order, the quadrupole 
formula predicts for the radial velocity f = — 16/5 (M/r)^ and for the ratio between the radial and tangential velocity 
r/{ujr) = -2/3cj/[j2 = -16/5 (M/r)^/^ = -16/5 {Mcu f/^ , where we use uj^r^/M = 1. In Fig. [S] we show how the 
above relations are satisfied by the numerical simulations. For simplicity we only consider the high-resolution run 
d = 19. Quite interestingly, the curves —16/5 (M/r)^/^ or — 16/5 (Mw)^/'^ average the behavior of r/{u;r) during 
the inspiral part and converge to it at later times. Between 30-50M before the formation of the CAH, we notice an 
abrupt change in the behavior of the ratio between the radial velocity and the tangential velocity, which suggests the 
presence of a blurred dynamical ISCO with subsequent plunge [1^. Even during the plunge, the radial velocity is 
still much smaller than the tangential velocity, reaching the value of 20% only at the end of the plunge. This result 
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Symbol 


Type 


Computed 


I^A 
WNQC 


orbital frequency 
-2C22 frequency 
-iCim dominant (circular-polarized) frequency 
-2C22 Newtonian quasi-circular frequency 


from AH centers (see Sec. IIB) 
by tracking wave peak (see Sec. IIB) 
from Eq. (l23ll 
from Eq. JSZl 



TABLE IV: We summarize several frequency variables used in the text. 



is a further confirmation that the numerical, equal-mass dynamics is quasi-circular until the end, as predicted by the 
EOB approach psj . 



IV. THE INSPIRAL 



The analysis in Sec. IIII Dl has shown that, despite the presence of an initial eccentricity, the dynamics is quasi- 
circular. If the dynamics is sufhciently quasi-circular, then it should be possible to model the inspiral waveform and 
frequency using Newtonian and PN methods together with the quadrupole formula. In this section, we will compare 
the numerical waveforms to the expected results from Newtonian theory and PN theory p6l | assuming an adiabatic 
inspiral. In a subsequent section (Sec. IVII[) we will also consider the non-adiabatic EOB model [27, 28, 29, 30] and 
Pade approximants [33| . Our analysis should be considered as a first-order attempt to assess the closeness of analytical 
and numerical results. More rigorous comparisons will be tackled in the future when numerical simulations start with 
initial conditions that more accurately model a binary on an adiabatic path closer to that describe by PN methods, 
and with simulations that have smaller or better understood systematic errors. 

In addition to examining the full waveforms, it is useful to focus attention on the angular frequency of the waves 
and the underlying orbital motion. From the evolved data, there are several methods for determining the orbital 
angular frequency. The most direct measure is obtained by tracking the coordinate locations of the centers of the 
AHs of each individual BH. We label this measure of the frequency by ujc- Because it is based directly on generalized 
harmonic coordinate values, this measure of lo is susceptible to gauge effects. A second method for determining lo is 
to track the phase of the maximum of ^"4 in the equatorial plane as it intersects the extraction surface at r = 50A/. 
We denote the orbital angular frequency determined by this method by lo\. Since the angular resolution at which 
\l/4 is sampled is coarser than the temporal resolution, we use spatial interpolation to find the phase (/>max(i) of the 
maximum at time then smooth the curve in t before computing u)\ = d(j)max/dt. As a third method, we note that if 
a complex signal f{t) has a dominant frequency and it is circular polarized, then that frequency is given by Im[///], 
where the dot (') denotes a time derivative. In terms of the mode amplitudes -2Cim{t,r), the dominant circular 
polarized frequency can be estimated by 



-Im 



2Clm 



(23) 



where we note that m in this equation is the azimuthal index and should not be confused with the total mass. These 
different definitions of the frequency are summarized in Table IIVI 

Figure [7| compares the orbital angular velocity M Ld{t) obtained by these approaches, for the initial separations 
d = 16 and 19. First, note that various frequencies have been appropriately shifted in time to account for the 
wave propagation time to the extraction sphere. The initial numerical waveform is dominated by spurious radiation 
associated with the initial-data, and lj\ and ijJd2 are quite noisy at early times. Though it is somewhat difficult to see 
in these plots, LjJ\ and ti;D2 also extend to earlier times than does Wc- This is a manifestation of the fact that uj\ and 
LJ-D2 are obtained from information at the extraction sphere at r = 50M. Also small and difficult to see, we note that 



Note that when comparing with analytical models we assume that the binary total mass M, introduced in Sec. Iml as the sum of the 
irreducible [6Cll | BH masses computed from the AH, coincides with the rest masses appearing in the PN waveforms, and is constant. 
In a numerical evolution, the mass estimated from the AH can change during evolution. In these simulations we believe most of this 
change is due to numerical error, though in principle part of it could be accretion of gravitational energy. Also, given that the AH is a 
coordinate dependent object part of the change could be gauge-related, though this is unlikely. Regardless of the source, for the highest 
resolution simulations the change in M is relatively small. For example, for the d = 19 case, we find that the maximum drift in M is 
less that about 0.2% before CAH formation. After the CAH, the final AH mass drops by about 1% in the last lOOM. Those variations 
are within or smaller than other errors present in the numerical simulation. 
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there are unexpected deviations at early times in ujc- We find that all three measures of oj agree quite well except 
near the beginning of the evolution and near the end of the inspiral 30-40 A/ before the time of the peak in \^4\. 

The aberrant behavior of uo at early times is primarily due to the use of conformally-flat initial data and, conse- 
quently, to a lack of physically realistic initial radiative modes. This strongly affects lox and ljd2- The small anomalous 
behavior in lJc at early times might also be caused by artifacts in the initial data, and note that imposing spacetime 
harmonic coordinates at the initial time does create some "artificial" coordinate dynamics. Each of these methods 
for measuring uj are based directly on generalized harmonic coordinate values, and are susceptible to gauge effects. In 
particular, the coordinate position of an AH is certainly not a gauge-invariant quantity, and given that the BHs are 
in the strong-field region of the spacetime there is no a priori reason to expect the coordinate locations to have any 
simple mapping to what one may describe as the physical orbit. It is therefore somewhat surprising how well these 
"almost-harmonic" coordinates describe the orbit — more examples of this are given in the next section. 



A. Newtonian quadrupole approximation 

Now consider a Newtonian binary in a circular orbit with orbital angular frequency lo. For a binary with reduced 
mass /i = mim2/M and mass ratio v = fi/M, the standard quadrupole formula yields 

rAf*4 32y^zy(A/u;)*^/3 g-2»M-0o)_2y22 + e2'('^*-'^n)_2i^2-2 (24) 

where 0o fixes the initial phase of the orbit and assuming right-handed rotation about the positive z-axis. If we 
replace ujt by the accumulated phase of the orbit 

= / ujit')dt', (25) 
Jo 

then we find that we can approximate £ = 2 modes of the inspiral waveform by 

_2C2^(t) = 32^1 1/ [Mcj(t)]8/3eT2^Wt)-0o), (26) 

Note that we have assumed an adiabatic inspiral and have replaced the constant orbital angular frequency uj of 
the circular orbit with a time-dependent orbital angular frequency uj{t). We refer to the result in Eq. ([26]) as the 
Newtonian quadrupole circular orbit (NQC) approximation. 

This result can be used in two ways. First, if we assume that Eq. (|26p provides a good approximation to the 
waveform, then we can extract uit) from the waveform via 

/ . I— \ 
Mu;{t)^i—^^\_,C,^{t)\\ . (27) 

In Fig. [7] we have also plotted wnqc obtained by the Newtonian quadrupole circular orbit approximation of Eq. (j27p . 
As noted previously, near the beginning of each evolution, artifacts from the initial data dominate the waveform and 
this leads to large inaccuracies in uj\ and wd2- Similar inaccuracies at early time are also seen in cjnqc- Near the 
end of the inspiral, when the orbital motion is no longer close to circular, we should not expect Eq. (|27p to yield an 
accurate value for AIuj and we see that the NQC method is systematically underestimating the value of Mui. 

The inconsistency of these methods for determining uj near the end of the inspiral should remind us that the 
very notion of "orbital angular frequency" becomes poorly defined after the formation of a common horizon, coc 
terminates near the end of the inspiral as a common horizon forms. The various methods agree quite well until about 
a quarter of an orbit before the formation of a common horizon (see Fig. |8] below). At this point, which we refer to 
as the "decoupling point", coc separates from cox and ujd2 and begins to rise more rapidly. Determining the precise 
point of decoupling is difficult due to numerical noise in the frequencies, but it seems to occur close to a value of 
Mojdec = 0.15 ± 0.01. Incidentally, this moment of decoupling also seems to coincide with the time the centers of the 
individual AH's cross what we estimate to be the late time co- rotating light ring of the final BH (see Sec lVI|) . Finally, 
in Sees. IIV Bl and IVIIi we will again examine the orbital angular frequency of the numerical models and find that co 
from 3PN- adiabatic and EOB circular orbits agrees well with ujc- 

The second way that the NQC approximation can be used is to estimate -2(^2^ (i) by using co extracted from the 
evolution. Figure [8] compares the real part of -2(^22(0 with the waveform estimated using the Newtonian quadrupole 
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FIG. 7: Orbital angular frequencies for the d = 16 and d = 19 cases evaluated using several different methods. The solid (black) 
line labeled uJc displays Muo as determined from the coordinate locations of the center of each BH's AH. The solid (green) line 
labeled uj\ displays Mlu as extracted by tracking the phase of the peak in 'I'4 at an extraction surface placed at r = 50M. The 
long-dash (red) line labeled uj-d2 displays the dominant frequency in -2C22 obtained using Eq. ((23} • The dash-dash-dot (blue) 
line labeled cjnqc displays the orbital angular velocity obtained from -2C22{t) using Eq. (|27l) . Finally, the vertical dotted (blue) 
line marks the approximate time that a common AH forms. 



circular orbit approximation of Eq. (pS)) for both the d = 16 and 19 cases. In the upper half of each plot, we use Uc 
for the orbital angular frequency. In the lower half of each plot, uj\ is used. We do not consider reconstructing the 
waveform from tJD2 or wnqc because these where themselves derived from -26*22 (i). 

A benefit of examining these plots is that they give a clear indication of how much of the initial waveform is 
contaminated by artifacts from the initial data. This can be most clearly seen in Fig. [5] in the comparison with lOc 
where we find that the estimated waveform begins later than the extracted waveform. The reason is that Wc is a 
function of "coordinate time" while the extracted waveform is a function of "retarded time" at the extraction radius. 
So, the beginning of the estimated waveform marks the earliest time that a waveform produced by the numerically 
evolved inspiral motion could begin. The numerical signal preceding this is due entirely to the unphysical initial 
radiative content of the initial data. This signal precedes the inspiral waveform because it originates from spatial 
locations in the domain that are closer than the binary to the extraction sphere. An initial segment of the true 
inspiral signal is also contaminated because of initial-data artifacts propagating to the extraction sphere from beyond 
the binary. If we make the reasonable assumption that the most significant contributions to the initial-data artifacts 
originate within the extraction sphere located at r = 50M, then we should expect a total of around 2 x 50M of the 
signal to be contaminated as measured in the retarded time of the extraction sphere. This number cannot be exact 
since we expect the initial-data artifacts to be strongest close to the center of the extraction sphere, and also because 
of variations in the coordinate speed of light in the strong field region. 

Both halves of the plots in Fig. [5] show a clear mismatch at early times. Because uj\ is constructed from information 
at the extraction sphere, it shows an initial pulse of radiation that is clearly an artifact of the initial data. However, 
the level of contamination of the wavform decays quickly following this initial pulse and appears to have become 
insignificant by a time of 30 A/ to 50M following this "initial-data pulse" ^. 

A striking feature of these figures is the excellent agreement between the estimated and extracted waveforms 
following the initial contamination and up to a short time before the formation of a common AH. During this phase 
of the inspiral. Fig. [1] clearly shows that the motion of the binary is not circular. Nor is it the smooth adiabatic 



^ We also note that if we were computing -2C22 without assuming the Keplerian relation LtP" jM = 1, the agreement would be better 
at earlier times because the Keplerian relation has its largest error there (see Fig. (5)1. 
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FIG. 8: Comparison of numerical and NQC inspiral waveforms for the d = 16 and 19 cases. In all plots, the solid (black) line 
displays -2C22{t) from the numerical waveform and the vertical long-dashed (blue) line marks the approximate time that a 
common AH forms. In the upper plots, the dashed (red) line displays -2C22(t) as computed from Eq. (|26|) using ujc and (j){t) is 
obtained from Eq. (|25p . In the lower plots, the dashed (red) line displays -2C22{t) as computed from Eq. (|26p using uli\. The 
longer evolution of the d = IQ run leads to significant changes in scale for -2C22{t). Therefore, we change scales for the final 
90 Af of the evolution. 



inspiral that we would expect from an astrophysical binary that has evolved from much larger separation. In fact, 
the observed motion exhibits a small radial oscillation about this "desired" motion. This effect is most easily seen 
in the longer d = 19 evolution, and also in the plots of oj in Fig. [7] (see also Fig. [51 and Fig. [TO] which includes a fit 
to the PN form for oj expected for adiabatic inspiral from large separation). The point we want to emphasize is that 
Eq. (j26p gives an excellent approximation for the waveform, even when the motion is clearly non-circular, so long as 
the phase (pit) and orbital angular velocity Lu{t) accurately incorporate the non-circular aspects of the orbital motion. 
As mentioned above, the NQC approximation appears to work quite well up until about 1/4 of an orbit (half a full 
wave cycle) before the appearance of a common horizon. 

Another rather intriguing example demonstrating the adequacy of the quadrupole formula, and how well adapted 
the numerical coordinate system is in describing the binary motion is shown in Fig. [5) For brevity we focus on a 
single example here, comparing the real part of the -2C22 component of the d = 19 waveform to the same component 
of a waveform calculated using the quadrupole formula ^ for two point sources of mass M/2 following trajectories 
given by the coordinate locations of the AH's from the simulation. The latter curve ends when a common AH forms, 
and has again been shifted in time by a constant amount to account for the propagation time for the wave to reach 
the extraction surface. The difference between this comparison and the preceding NQC comparison is we have not 
assumed circular orbits, using instead the detailed orbit motion obtained from the simulation. Not only does the good 
agreement testify to the well-suited nature of the coordinates, it shows that the quadrupole formula does a remarkably 
good job of capturing the dominant physics of GW emission during the entire merger regime prior to common AH 
formation. 



Here we mean taking directly four time derivatives of the coordinate motion of the centers of the individual AHs. 
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FIG. 9: Comparison of the real part of the -2C22 component of the waveform from the d = 19 (ft/2) numerical simulation 
(solid curve) versus the same quantity calculated by using the quadrupole formula (dashed curve) for two point particles, each 
of mass M/2 and following trajectories given by the coordinate motion of the AH's in the simulation (see Fig.[T|. The waveform 
from the orbit ends once a common AH forms. 



B. Adiabatic post-Newtonian model 



The post-Newtonian approximation to the two-body dynamics of compact objects provides the most accurate 
predictions for the motion and the GW emission during the inspiral phase, when the weak-field and slow-motion 
assumptions hold. 

In a more rigorous analysis we would compare the numerical and analytical dynamical quantities expressed in the 
same coordinate system and gauge. Here and in Appendix [B] we limit the comparison to a few gauge-invariant 
dynamical quantities, such as the orbital frequency the orbital phase ^, the energy flux and the angular- 
momentum flux i^j, expressed in terms of the instantaneous orbital frequency and/or the time of a stationary observer 
at infinity. This is also motivated by the fact that previous studies [2^, [33| have shown that PN-approximants to 
dynamical quantities are more robust (under change of PN order) if expressed in terms of gauge-invariant quantities, 
notably the instantaneous orbital-frequency Mlo. At the present time, PN calculations provide the orbital frequency 
through 3.5PN order [2^ if spins are neglected, and through 2.5PN order [t^I if spins are included. As mentioned 
in Sec. [TT] and shown in Table HI the numerical initial data describe BHs which carry a small spin aligned with the 
direction of the orbital angular-momentum £. For this reason we include spin effects in the PN approximants. 

In this section we limit the analysis to the so-called PN-adiabatic model. In Sec. IVIII we shall investigate the 
comparisons with the EOB model which goes beyond the adiabatic approximation. In the PN-adiabatic model the 
waveforms are computed assuming that the motion proceeds along an adiabatic sequence of quasi-circular orbits. More 
specifically, one assumes f = and evaluates the variation in time of the orbital frequency to from the energy-balance 
equation dE/dt — —Fe, where E is the two-body energy and Fe is the GW energy flux. In particular, E and Fe are 
first computed for circular orbits and written as a power expansion in AIoj, then uj{t) — ~Fe{uj)/ {dE{u})/doj). The 
adiabatic sequence of circular orbits ends at the conservative Innermost Circular Orbit (ICO), i.e., the ICO evaluated 
from the conservative dynamics by imposing (dE / duj)ico = [1^. The study of Ref. (see in particular Figs. 4 
and 5 and discussion around them) and Ref. [53|, showed that waveforms computed in the adiabatic approximation 
(which are very accurate at large separations) can have a non-negligible phase difference with respect to waveforms 
computed in the non-adiabatic approximation, even before reaching the last stable orbit. The accuracy of our numer- 
ical simulations and the nature of the initial data will not allow us to explore these phase differences. Nevertheless, 
we have found it useful to use the adiabatic PN model as a diagnostic of the last few cycles of the numerical evolution. 

As discussed in Sec. IIV Al it takes a certain time for the evolution to settle to a quasi-circular orbit. Moreover, the 
numerical results contain a non-negligible amount of eccentricity. For these reasons we shall evaluate the PN-adiabatic 
approximant which best averages the numerical orbital frequency until either the dynamical ISCO, the decoupling 
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FIG. 10: We compare the NR and three analytical orbital frequencies obtained by fitting (i) both tc and x (dashed line), (ii) 
only tc and using the nominal x value from Table |T] (dot-dashed line) and (iii) only tc and assuming x = (dot line). The left 
panel refers to the run d = 16, the right panel to d = 19. From top to bottom, the continuous horizontal lines mark the ICO 
frequencies as predicted by the 3PN-adiabatic energy and the decoupling frequency. The dashed horizontal line in the right 
panel span the frequency range when a dynamical ISCO could be present. 



frequency, or the CAH. Again, these issues will be overcome when numerical simulations starting at larger separation, 
and from initial conditions that more accurately model an adiabatic inspiral, become available. We notice that in 
principle there could be non- negligible differences between the instantaneous orbital frequency as defined in PN theory 
and in the numerical simulation. In the latter, the orbital angular frequency is calculated from the coordinate locations 
of the centers of the AHs of each individual BH. However, since in Sec. lIVIBj we have found that the numerical orbital- 
frequency agrees quite well with the numerical GW frequency extracted at larger radii, we expect that the differences 
are small. 

Defining t = v [tc — t)/5M, where tc is the time at which the orbital-frequency diverges (time of coalescence, not 
to be confused with the decoupling time) we have 
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where 5M = mi — TO2, the spin variables are 



S = Si + S2, 



(29) 
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FIG. 11: We compare the NR and analytical Re[-2C22] obtained by fitting the orbital frequencies. The fit is done for (i) both 
tc and X (dashed line), (ii) only for tc and using the nominal value for x from Table |I] (dot-dashed line) and (iii) only for tc and 
assuming x = (dot line). The left panel refers to the run d = 16, the right panel to d = 19. The dashed vertical lines span 
the region during which a dynamical ISCO could be present. The waveform cycle beyond this point could be associated with 
the plunge cycle. 
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with S,; = x^m^Si. In Eq. 
orbital angular-momentum £ 



we have denoted with Si and the spin components along the direction of the 
The orbital phase through 3.5PN order reads 



1 



T-5/8 



/3715 55 



I' { V8064 
47 Sg 75 SM 
16 M2 ^ 6417 M2 
9275495 284875 



96 

^1/4 

1855 



T-3/8 



rl/4 



14450688 258048 
/ 508265 3665 
V 258048 7168'^ 
831032450749357 _ 
57682522275840 ~ io" 
126510089885 2255 



rl/8 



2048 
Si / 37265 

53 



38645 65 , , , , 

1 ly TT In r 

172032 2048 ' ^ ' 



V 57344 
107 



1735 
7168' 



SM 



ln(r) 



56 



-C- 



— 1 (—] 
448 "^1256/ 



4161798144 



2048 



154565 
1835008^ 



1179625 
1769472' 



-1/8 



188516689 488825 



173408256 516096 



141769 
516096' 



1.2 ] 7rr"i/4 



(31) 



where C = 0.577 • • • is the Euler constant and (t>p is an arbitrary constant. The non-spin terms in Eqs. (|28p and 
([3T|) are given by Eqs. (12) and (13) of Refs. [1^, [7l|, while we evaluated the spin terms through 2.5PN order using 
the recent results of Ref. [t^ (we use the constant spin variables as defined in Ref. [70|), and we neglected spin-spin 
contributions. 

For each of the three runs, we determine the time of coalescence tj, and the spin-magnitude X = Xi = X2 by fitting 
the PN orbital frequency (j28p to the numerical orbital frequency using a non-linear least-squares method. Figure fTUl 
shows the results for d = 16 and d = 19. Due to the initial burst of radiation related to the initial conditions, we 
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FIG. 12: We plot the numerical h+ (left panel) and hx (right panel) extracted along the direction perpendicular to the orbital 
plane (dark continuous line) and compare them with PN-adiabatic predictions with the phase evaluated at 3PN order and the 
amplitude evaluated either at Newtonian order (dashed line) or through 2PN order (light continuous line). These results refer 
to the d = 19 run. 



remove from the numerical data the first « 22M and « 34M, for the d = 16 and d ~ 19 runs, respectively. In 
Sec. IIIIDI we discussed the possible presence of a blurred dynamical ISCO [11] occurring 30-50M before the CAH 
forms, at Mwdynisco = 0.078-0.097. The orbital frequency at the conservative ICO evaluated at 3PN order with 
X = 0.080 is M wico = 0.143 and with x = 0.063 is M wjco = 0.140. We mark all these frequencies in Fig. [TU] along 
with the decoupling frequency. In principle, the adiabatic PN waveform should be used only until the last stable 
orbit, since it was derived from the balance equation for circular orbits which ends at this last stable orbit. However, 
since the two-body motion predicted by the numerical simulations is rather adiabatic and quasi-circular until the 
CAH forms [28^, we extend the PN waveforms through the plunge until almost that point. More specifically, we fit 
until the time at which the numerical orbital frequency decouples from the GW frequency, Mcudcc ~ 0.15 ± 0.01. 
Notice that by fitting the time of coalescence tc and x '^6 ^-re fitting the initial value of the orbital frequency. We 
find that the 3PN-approximant best fits the data with x'*^^^ = 0.111 and 1"^=^^ = 220.1 M [M uj{t = 0) = 0.04789], 
;^^d=i9 ^ 0.0874 and 4=^^ = 509.1 M [M u;{t = 0) = 0.03477]. If we fit until the CAH time, we find x'^^^^ = 0.0812 
and 4=^^ = 217.6 M [M iv(t = 0) = 0.04715], x'^^^^ = 0.0626 and t'^=^^ = 506.2 m [M u}{t = 0) = 0.03447]. Those 
values are closer to the nominal x values of Table HI However, we think this is accidental. Finally, notice that if we fit 
until the dynamcal ISCO Mwdynisco = 0.096 we find x'^^^^ = 0.113 and t^=^'^ = 514.3 M [M iv(t = 0) = 0.03502]. 

In Fig. [TO] we also show curves evaluated using the nominal x value of Table H] To understand how spins affect the 
PN-adiabatic orbital frequency, we show in Fig. [TO] also the case in which we fix x = and fit only tc- The latter 
values produce a difference of ~ 0.3 GW cycles at the CAH time with respect to the case where we fit both tc and x- 
Although we can use the numerical results to discriminate between several PN-adiabatic models with spin, it is not 
clear which role the spin variable is playing in fitting the data. In fact, the spin values obtained from the fit are also 
affected by the eccentricity present in the numerical data but absent in the analytical model. 

Using the time of coalescence and spin values obtained from the fits, we plot in Fig. [TT]the waveforms ReJ_2C'22]. 
Notice that the GW phase differences between the fitted models is smaller than the maximum GW phase error 
estimated in Table iHll using lower resolution runs. This gives a concrete example of how cumulative phase error in 
a numerical simulation translates to uncertainties associating PN model parameters with the numerical waveform, 
despite the deceptively small phase error after maquillage. PN-adiabatic models can also fit the data of lower resolution 
runs and give initial orbital frequencies larger than those found for the high resolution runs. 

So far, when comparing with numerical waveforms, we have neglected higher-order PN corrections to the GW 
amplitude and have restricted the comparison to Re[_2C22], i.e., we used waveforms in the so-called restricted ap- 
proximation. In Ref. 't?] the authors evaluated the ready-to-use PN waveforms /i+ and hy through 2PN order in the 
amplitude (see Ref. jTSj] where this computation has been pushed through 2.5PN order). In Figs. [TO] we compare the 
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FIG. 13: In the left panel we plot the orbital frequency for several PN approximants computed from Eqs. (|28p and (|32p whose 
frequencies coincide at time t = 0. In the right panel we plot the differences between the number of GW cycles at 3PN computed 
from the analytical Eq. (I28|l . and at nPN order computed from Eq. (|32p . The initial frequency coincides with the one computed 
from the fit of the d — 19 run when x — 0. 



numerical h+ and hx , which are obtained by integrating -04 twice in time, with the analytical /i+ and hx as given by 
Eqs. (2)-(4) of Ref. [72] ■ The phase is computed from Eq. ([31]) at 3PN order with the values of and x obtained from 
the best fit. The waves are extracted along the direction perpendicular to the orbital plane. Because the BH masses 
are equal, only the 2"'^ harmonic is present. We would conclude that higher-order PN amplitude corrections have a 
mild effect in the waveform emitted by equal-mass binaries. However, we notice an oscillating behavior in the PN 
approximation to /i+ and hx- For example the IPN correction is rather large and opposite in sign to the Newtonian 
correction, resulting in a significant reduction of the amplitude of the signal. The 1.5PN and 2PN corrections undo 
this effect. This oscillating behavior seems to also affect the higher multipoles Cim- In fact, we checked that C44 is 
well approximated by the 3PN-adiabatic model for the phase, if computed at (leading) IPN order in the amplitude, 
but the agreement become worse when 2PN corrections are added. We plan to investigate in more detail the effect of 
higher-order PN corrections and higher multipoles in the future. 

To obtain more robust comparisons between PN and numerical predictions, it would be preferable to start the 
numerical evolution where we are confident that the PN expansion can be safely applied. First, we notice that if 
we were computing w as a function of time instead of from Eq. (I28p . then integrating numerically the following 
equation [2£ 
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which can be derived expanding Co — ^ F [lo) / {dE / cLo) in powers of {Mlj), we find some differences from Eq. (j28p . 
In particular, at 3.5PN and 2.5PN orders, lo derived from Eq. ([28|) reaches a maximum and then starts decreasing, 
becoming negative. By contrast, this behavior does not occur in the u derived by numerically integrating Eq. (j32p . 
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FIG. 14: We plot the differences between the number of GW cycles at 3.5PN and at nPN order versus the 3.5PN uj. All quantities 
are computed integrating numerically Eq. (|32|l with spins set to zero. The initial and final frequencies are Mu) = 0.004 and 
Muj — 0.026 in the left panel, Mui = 0.02 and Mcj = 0.13 in the right panel. They correspond to an equal-mass binary 
sweeping in the most sensitive frequency band of LIGO from ~ 43 Hz to ~ 280 Hz of mass (3 + 3)Mq in the left panel and of 
mass (15 + 15) in the right panel. 



In Fig. 1131 we show the differences in the orbital frequency and the number of GW cycles if the latter quantities were 
computed from Eq. ([5^ at different PN orders but with the same initial frequency luq — 0.03361. We consider here a 
non-spinning binary. In the right panel of Fig. [13] we compute the differences in the number of GW cycles between the 
3PN UJ from the analytical Eq. and several PN uj from Eq. ([5^ . Quite interestingly the 3.5PN order computed 
numerically is very close to the 3PN order computed analytically, whereas the 3.5PN order computed analytically has 
almost half a cycle of difference at the end of the inspiral. These differences are a consequence of the fact that at 
such (close) initial separation, M uj ^ 0.033, the differences between PN-approximants are not negligible. This fact 
is better illustrated in Fig. [TH We start the evolutions at Mw = 0.004 (left panel) and M uj = 0.02 (right panel) 
and plot the differences between the number of GW cycles at 3.5PN and at nPN order versus the uj computed at 
3.5PN order. All quantities are obtained by numerically integrating Eq. (|28p with spins set to zero. 2PN and 2.5PN 
approximants accumulate large differences from the 3.5PN-approximant when evolving from larger separations (left 
panel) . 

Thus, summarizing, the phase computed at 3PN order from Eq. or at 3.5PN order from Eq. ([5^ best-fit the 
numerical results. If we wanted to investigate to which PN order, notably 3PN or 3.5PN, the NR orbital frequency 
is closest to, we would need to start the numerical evolution at frequencies smaller then the one used in the o? = 19 
run, which is 0.033. 



V. THE RING-DOWN PHASE 

During the ring-down, the GW can be decomposed in terms of the quasi- normal modes (QNM) of Kerr[67i[7l.[75t. 
These modes are distinguished by their longitudinal and azimuthal indices i and to, as well as by their overtone 
number n. Each mode has a particular frequency uJimn and decay constant r^mn which are functions of the Kerr 
parameter a and total mass Mf of the background BH that is being perturbed. To shorten the notation, we will 
introduce the complex frequency toimn and use * to denote complex conjugation: 



Wfmn = '^Imn — i/Timn- 



(33) 
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Following Ref. [7g|, the ring down can be expressed, in terms of the Weyl scalar as 

where Cemn, C[^^, 4>imn, and are real constants, and Stmn = Sira{ailiimn) are the spin-weight -2 spheroidal 

harmonics^ implicitly evaluated at the complex QNM frequencies. The primed terms are necessary because, for a 
given (€, m, n) and a fixed non-vanishing angular momentum, there are two solutions of the eigenvalue problem. To 
make the notation as clear as possible, we will always take the real frequencies and decay constants to be non-negative, 
so Wfmn ^ and r^mn ^ 0- Of the two solutions to the eigenvalue problem for fixed {^, m, n), one solution has positive 
frequency and one negative, and the complex frequencies are related by —u>imn — '^i^m (^^^ -R-^f. i76j] for a full 
discussion). Because of this relationship, it is only necessary to determine the positive (or negative) frequency modes. 
In Ref. [76| the authors compute the positive frequency modes and choose the convention that ujimn ^ ^i-^nm with 
equality in the case that m = or a = 0. Finally, we note that with these conventions, it is necessary to introduce an 
overall sign change on the real frequency in the equations of Ref. (76| in order for the signs of the frequencies of the 
various modes to agree with numerical simulations. 

The decomposition of \E'4 in terms of spin- weight -2 spherical harmonics is given by Eq. (jl7p . In order to relate 
the expansion coefficients in Eqs ([M)) and fTT]) , we need the expansion of the spheriodal harmonics in terms of the 
spherical harmonics. Following Press and Teukolsky [TSj . 

Stmn = ■Au"mn-2Yl"m- (35) 

Using the orthonormality of spin-weighted spherical harmonics, we find that 

l"n 

So, in principle, a spherical harmonic mode amplitude -2Cem of the ring-down signal can contain a contribution from 
any of the negative frequency modes with azimuthal index m and from any of the positive frequency modes with 
azimuthal index — m. Note that in the second version of this expansion, the complex coefficients Awmn have been 
absorbed into the new real expansion coefficients Cu"mm C'ge"mn^ 4>U"mn, and 4''u"mn^ where each coefficient now has 
four indices. 

The expansion coefhcients Awmn depend on the product of the Kerr parameter and the complex QNM frequency 
aujtrnn and for sufficiently small values can be determined via perturbation theory (cf. Ref. [ill). For example, using 
first order perturbation theory, we find 

-2^22 - ^{c22„e-*('^^^"*+*^^") +C^_2ne'^"=*-^"*+'^^-^"^ (38) 

n 



+ ^a(i32„(4 + acD32„)C32ne-'('^-"*+^^-) (39) 

- ^\[\au^UJA + ac2;*_2„)C^_2„e^'^--*+*--' (40) 
+ ^(aci42„)'C42„e-'(^-"*+^-") (41) 

+ ^(«'^^2J'C^2„e^'--*+^--)}. (42) 

We have extracted the various QNM contributions to the _2C22(t) ring-down signal in the following way ®. At late 
times, we expect the i! = 2,m = 2, n = QNM to dominate. We fit the signal after time tr + ipcak to this single mode 



We note that the 9 depencence of spheroidal harmonics is connected to the separabihty of the Kerr metric in terms of Boyer-Lindquist 
coordinates. While our spherical coordinate system is not Boyer-Lindquist, the differences are not significant in the wave zone where 
the waveform is extracted. 

* After the work to fit the ring-down modes was completed, a similar approach was posted in the preprint archives |77l| . 



23 



TABLE V: Richardson extrapolated values for the angular momentum ratio a/A4f and final mass ratio Mf/M from the ring- 
down fits to Re[_2C22] for the d = H, d = 16, and d = 19 cases. For each separation we provide estimates when A'' overtones 
of the f = 2, m = 2 QNMs are used. 





d = 


13 


d = 


16 


d = 


19 


N 


a/Ms 


Mf/M 


a/Mf 


Mf/M 


a/Mf 


Mf/M 





0.724 ± 0.002 


0.948 ± 0.006 


0.72 ±0.01 


0.945 ± 0.003 


0.702 ± 0.007 


0.946 ± 0.002 


1 


0.723 ± 0.002 


0.945 ± 0.008 


0.72 ±0.01 


0.942 ± 0.001 


0.706 ± 0.004 


0.947 ± 0.001 


2 


0.735 ± 0.002 


0.955 ± 0.006 


0.725 ± 0.007 


0.946 ± 0.002 


0.711 ±0.002 


0.949 ± 0.002 


3 


0.732 ± 0.006 


0.951 ±0.003 


0.725 ± 0.007 


0.946 ± 0.002 


0.709 ± 0.003 


0.947 ± 0.001 



using non-linear regression and choose tr to minimize the error in the fit. There are four dimensionless parameters in 
this non-linear fit: C22O) '/'2201 "11^220? ^^^id T^inl'tn. However, instead of fitting directly for these four parameters, we 
treat mujimn and Timn/'Tn as functions of a/Mf and Mj /m which can be obtained via interpolation from tabulated 
values (cf. Tables II-IV of Ref. [76]) or via approximating functions (cf. Tables VIII-X of Ref. [zl])- The advantage 
of using {a/Mf, Mf /m, €220, ^220) for the set of fitting parameters comes when we fit to additional modes. 

If we knew a/Mf and Mf/M precisely from the fit to the dominant mode, then we could directly compute the 
values for Mujgjrm and Temn/M for all additional contributing modes. All that would be necessary for determining the 
contribution of each additional mode would be to fit for its amplitude and phase. However, we will not know a/Mf 
and Alf /M precisely from the fit to the dominant mode. Therefore, we treat the frequency and damping constant for 
each mode as functions of a/Mf and Mf/M so that they are determined consistently when we fit multiple modes to 
a given signal. 

To make our procedure more explicit, we note that for the -2C22(^) ring-down signal, the contributions of modes 
with £ > 2 seem to be near the level of numerical precision, as are the positive frequency modes with £ = 2 and 
m = —2. Therefore, we take as our fitting function: 

JV 

-2C22 = ^C222rie"*'^'^^^"^"/*^^'*^^/™^^*"*P'="''^+'^^^^"\ (43) 
n=0 

We fit separately the real Re[_2C22] and imaginary Im[_2C22] parts of the -2(^22 (i) ring-down signal without any 
phase shifting of the numerical waveform. Separate fits were performed because simple non-linear least-squares fitting 
was used. During each fit, numerical data for the waveforms for times prior to + tpcak are removed from the signal. 
Starting with TV = 0, we fit the data for a sequence of values for tr and choose as our final the value that produces 
the smallest error estimate for a/Mf and Mf/M. We include additional overtones {N > 0) successively, using results 
from iV = fits as seeds for the = 1 fits, and so forth. For each value of N, we refit the entire function, so for 
A^ = there are 4 parameters in the fit, for iV = 1 there are 6, for A^ = 2 there are 8, and so forth (see Tables rvTIl 
IVIIIl and IIXI for and explicit list of the parameters being fit for each TV). At each set, we determine new values for 
a/Mf and Mf/M that are used consistently for all of the modes. Also, each time we include a new mode in the fit, 
we also fit the data for a sequence of values for tr and again choose our final tr for that set of modes by the value 
that produces the smallest error estimates for a/Mf and Mf/M. 

Tables IVlIl IVHIl and IIXI in Appendix |D] display the fit parameters for the -2C22{t) waveform obtained from initial 
data with separation parameter d = 13, c? = 16 and d = 19. These tables show the fits for only the high resolution {\h) 
runs and give result to 3 significant figures. In most cases, the errors in the fit suggest that only two significant figures 
can be trusted, but we display the additional digit in order to clearly illustrate the level of consistency in the fits. 
However, even if the accuracy of the individual fits were higher, we must still take into account the discretization error 
when estimating the value of parameters from the fits. Using Richardson techniques, we can for example estimate the 
value and error of the angular momentum and mass of the BH at the end of the ring-down phase by using fits to the 
medium (|/i) and high (^/i) resolution runs. Table IVl shows the results of this analysis for the angular momentum 
ratio a/Mf and final mass ratio Mf/m and includes results for the d = 13, d = 16, and d = 19 cases. There is 
considerable consistency in the value of the final mass ratio, with Mf/M « 0.95 for all separations. However, there 
is a discernible decrease in a/Mf as the separation increases. In fact each case differs by about 0.01 in value from its 
neighboring separation. This variation in the final spin of the coalesced BH is, in fact, completely consistent with the 
change in the spins of the initial corotating BHs. If we compute the total angular momentum contained in the spin of 
the individual BHs S, then we find respectively for the d = (13, 16, 19) cases S/M'j = (0.06, 0.04, 0.03). 

Figure [15] shows the quality of the fit to Ke[^2C22{t)] for the cases A^ = 0,1,2,3 and for the separations d = 16 
and 19. We note that by including modes through the n = 3 overtone, we can fit the ring-down quite well to times 
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FIG. 15: Comparison of numerical and QNM -2C22 ring-down waveforms for d = 16 and 19. All plots show the numerical 
ring-down waveform as a thin solid (black) line. The thick solid (red) line displays the fit of the ring-down signal using the 
first N £ — 2,m — 2 overtones beyond the fundamental. For plots containing A*' > overtones, we also include the fit residual 
from the previous value of A'^. This is displayed as a dashed (blue) line. The coefficients for the displayed fits are found in 
Tables IVlTIl and HXl 



preceding the point where |_2C22| reaches its peak. For each case beyond the fit to the fundamental n — mode, we 
include the residual of the previous fit. To be explicit, the residual displayed for iV = 1 is defined as the difference 
between the numerical signal and the fit obtained using the fundamental mode. The residual displayed for N = 2 
is the difference between the numerical signal and the n = 0, 1 modes used in the N — 1 fit. This residual gives an 
estimate of the remaining signal that is being fit. However, it is important to remember that for each value of N, 
the entire signal is actually being fit, including a redetermination of a/Mf and Mf /M for all the modes. The most 
important point to notice from the residuals is that for each value of A'' there is a clear signal that is being fit. 

A close examination of Tables IVIIHIXI reveals a significant level of consistency to the fits. For each separation d, 
the spin and mass ratios remain very consistent and the C22n and 4>22n coefficients remain quite consistent, as we 
increase the number of overtones included in the fits. This is true individually within the separate fits of the real and 
imaginary parts of -2C22, and consistency is also seen between the fits of the real and imaginary parts. While the 
t — 2, 171 — 2 QNMs seem to dominate the ring down signal in -2C22, the £ — 2, m = —2 modes and the modes with 
£ > 2 should be present. However, the remaining residual after the iV = 3 fit (not shown in any figure) has very low 
amplitude at times after the peak in |_2C22|. While there are some hints to structure, there is insufficient signal and 
the simple approach we have used for fitting does not yield consistent fits when additional modes are included. 

However, if we fix the values for a/Mf and Mf/M to the values obtained from the £ = 2, m = 2 fits, we can fit 
for the Cg±zn and 4)g±zn coefficients for a range of modes. Doing so, we find that the fundamental QNM with £ = 3, 
m = —2 has the most significant contribution, followed by the i — A, m — —2 and i — m — 2 fundamental modes 
at roughly comparable levels. Unlike the case of fitting only the £ — 2, m = 2 modes, adding in higher overtones 
when an increased spectrum of modes was considered did not lead to consistent fits. Part of the difficulty in finding 
consistent fits to the subdominant modes is likely due to the fact that the signal associated with these modes is close 
to the level of numerical precision in the waveform. However, it is also likely that more sophisticated fitting methods 
are needed. In particular, it would be useful to fit the real and imaginary parts of the waveform simultaneously. It 
may also be helpful to fit several -2Cem modes simultaneously. 

While fitting multiple modes is problematic in some cases, it is essential in others. For the case of -2(^32 (t), the 
dominant QNMs include both 1—2 and ^ = 3, both with m — 2. In fact, it was not possible to fit the ring-down 
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FIG. 16: Comparison of numerical and QNM -2C32 ring-down waveforms for d = 16 and 19. All plots show the numerical 
ring-down waveform as a thin solid (black) line. The thick solid (red) line displays the fit of the ring-down signal using the first 
N £ = 2,m = 2 and £ = 3,m = 2 overtones beyond the fundamental. For plots containing A'' > overtones, we also include 
the fit residual from the previoius value of N. This is displayed as a dashed (blue) line. The coefficients for the displayed fits 
are found in Tables IXj and IXU 



signal of ~2C32(t) without fitting simultaneously for these two modes. To be explicit, we take as our fitting function: 



N 



-2C32 = ^|C322„e"'('^"""(''/*^^'^'^^/^'^)(*"*p'='"''+'^^22"^ 



(44) 



71=0 



Fitting proceeds as with _2C22, starting with = and then adding successive overtones which allow us to fit to 
successively earlier times in the ring down. 

Tables 1x1 IXIl and IXIII in Appendix ID] display the fit parameters for the -2C'32(t) waveform obtained from initial 
data with separation parameter d = 13, d = 16 and d — 19. Figure \W\ shows the quality of the fit to Rc[_2C32(i)] 
for the cases A^ = 0, 1, 2 and for the separations d — 16 and 19. We note that by including modes through the n — 2 
overtone, we can again fit the ring down quite well to times preceding the point where I-2C22I reaches its peak. As 
for Ke[^2C22{t)], we also include the residuals of the previous fit. We note that the level of consistency of the fits, 
though significant, is not as high for -2^32 as seen for _2C22. 

During the ring-down phase, it is possible for a few percent of the final mass Mf and angular momentum aMf to 
be radiated away from the system. The Kerr QNM frequencies and decay constants are computed assuming that the 
mass and angular momentum they carry away constitute a negligible perturbation on the system. This raises the 
question as to whether or not the radiated energy and angular momentum are affecting the QNM fits. This issue will, 
of course, become more significant as the fits are pushed to earlier times. As we have seen for the cases of -26*22 and 
-2^32, fitting to earlier times in the ring down requires the use of higher overtones (n > 0) with shorter decay times. 
Because these higher overtones dominate the waveform only at earlier times in the ring down, we should expect some 
increase in the level of uncertainty in the fits as we incorporate these overtones. 

Finally we want to revisit the plots of the orbital angular frequency displayed in Fig. [7l The ux and uj]j2 frequencies 
continue beyond the inspiral phase and through the ring down. Beyond the inspiral phase, this frequency clearly cannot 
be associated with the orbital angular frequency. Rather they are half of the dominant GW frequencies seen in -26*22. 
In Fig. [T71 we plot this dominant frequency from a time about lOM before the formation of a common AH and 
through the ring down. In this range of times, lOc clearly decouples from ux and ujd2. As the dynamics transitions 
from the inspiral phase, the dominant frequency rises very rapidly, finally reaching a plateau associated with the 
dominant QNM ring-down frequency. Both ujx and lou2 agree quite well through both the transition and ring down, 
but we note that uj]j2 shows an unusual "beating" of the frequency during the ring down. 



26 




FIG. 17: Dominant frequencies during the ring down for the d — 16 and 19 cases evaluated using several methods. The short 
dashed (red) line is ujc- The dot-dashed (blue) line is ujx- And the thin solid (purple) line is ujd2- These three lines were 
previously displayed in Fig. [T] The thick solid (black) line shows the dominant frequency of Eq. (143 |l as fit to the real part of 
the -2C22 ring-down signal. The long-dashed (magenta) line shows the dominant frequency of a similar fit, but consistently 
using the both the real and imaginary fits to -26*22. The vertical dotted (blue) line marks the approximate time that a common 
AH forms. 



We also plot in Fig. [T7]the dominant frequency as measured by the fits to the ring down. Using the ring-down fit 
function in Eq. together with the fit value given in Tables IVntilXI yields an analytic expression for -2C22 through 
the ring-down phase. Because we independently fit Eq. (|43p to the real and imaginary parts of -2C22{t) we have three 
different ways that we can construct _2C22. We can take the coefficients for the fit from either Re[_2C'22] or Im[_2C'22] 
and use that set of coefficients exclusively in Eq. (|^ . Using the analytic representation of -2C22{t) we can compute 
the dominant frequency using Eq. (|23p with m = 2. A plot of this frequency using the coefficient obtained form the 
fit of Rc[„2C22] is shown if Fig. [17] with the label "Re RD fit (N=3) ." Notice that the frequencies agree well during 
the ring-down phase and show a period following the peak in I-2C22I where the frequency increases before reaching 
its plateau. Also, we see no evidence of the beating seen in 1^02 • 

However, if we instead construct an analytic representation of -2*^22 using the coefficients from the fits to both 
Re[_2C22] and Im[_2C22], we recover the beating of the frequency. To be clear, using Eq. (|23|) to construct the 
dominant frequency incorporates both the real and imaginary parts of _2C22. If we consistently use the coefficients 
from the fit to Re[_2C22] when constructing an analytic representation for Re[_2C22] and use the coefficients from 
Im[_2C22] for its representation, then we obtain the line labeled "Re+Im fits" in Fig.[T71 The plot of this frequency 
clearly shows a beating of the frequency and this is due to a small mismatch between the real and imaginary fit 
coefficients. It seems clear that the beating we observe in wd2 is caused by a similar effect. Essentially, numerical 
error is leading to a non-physical mode that is not circularly polarized, leading to the mismatch seen between the real 
and imaginary parts of 5*4. 



VI. THE (PLUNGE AND) MERGER 



In Sec. nil Dl we discussed the possible presence of a rather blurred dynamical ISCO which marks the beginning 
of the plunge phase. The latter ends when the CAH forms. The plunge has a duration of 30-50M, corresponding 
to 1-1.5 GW cycles. The plunge cycle has a slightly different shape than the inspiralling cycles when viewed in 
Re[_2C'22] (see Fig.[TT]), but it can barely be distinguished from the inspiraUing cycles when viewed in /i+ and (see 
Fig. [T2| . Quite interestingly we notice that the onset of the plunge phase seems to happen soon after the "knee" in 
the frequency curve (see Fig. [7|), and when the first change in the slope of the GW energy flux occurs (see Fig. 
especially the right panel) . A second change of slope in the frequency and GW energy flux seems to happen roughly 
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around the CAH, the third change occurs at the peak of the radiation. 

In Fig. [T5]we illustrate some other interesting features of the inspiral to ring-down transition, i.e., the binary BH 
merger. We plot the frequencies Uc and u}\, and the GW energy flux (multiplied by 100). Circles mark the position 
(time and frequency) at which the CAH forms and show when the coordinate separation between the BHs become less 
that the estimated co-rotating light ring of the final BH. The latter are coordinate dependent quantities. The light 
ring is an unstable circular null geodesic of the Kerr geometry in the equatorial plane of the BH, and we estimate the 
position of the light ring by noting that for a Kerr BH with a = 0.70, in Boyer-Lindquist coordinates the co-rotating 
light ring is a radial distance of w 1.17 times that of the outer horizon |39j. For an estimate of the light-ring location 
in the generalized harmonic coordinates of the simulation we took the late-time coordinate radius of the final AH 
multiplied by 1.17. Of course the different coordinates used to arrive at this value make it a rather rough estimate, 
though given that the notion of a light ring is not well-defined near coalescence it would not add much if we found 
the exact location. When the equal-mass binary reaches the CAH, only half of the total energy has been released. A 
little while before this point we observe that u\ decouples from Wc- Note also that at the peak of the radiation, 23% of 
the total energy has yet to be released. Here and in the following, the total energy refers to the energy radiated from 
the beginning of the simulation until the end. Thus, it doesn't include the energy radiated during the long inspiral 
preceding the initial time of the simulation. 

The results obtained in Sec. |Vl in particular Figs. [17] and discussion around them, suggests that the GW emission 
soon after the peak of radiation is caused by the excitation of the QNMs of the final Kerr BH. The frequency of the 
least damped QNM / = 2, m = 2, n = is responsible for the plateau, and the higher overtones (Z = 2, m = 2, n > 0), 
which have smaller frequencies and smaller decay times, should be responsible for raising the frequency from the peak 
of the radiation to the plateau. It is still not completely clear to us whether the higher overtones and/or other QNMs, 
e.g., I = 2,m 2,n > 0, are able to smoothly connect the decoupling frequency with the plateau. In fact, soon after 
the decoupling, there could be a very short non-linear phase, perhaps with strong mode mixing, that would preclude 
a description in terms of QNMs. To clarify those issues, it would be interesting to compute the binary BH metric 
around the decoupling point or the CAH, decompose it as a single Kerr metric plus perturbations, e.g., as done in the 
close limit approximation (tqI. [soI. [Sll. [s^ . and determine more precisely when the perturbative regime starts. In the 
next section, following a more phenomenological approach aimed at providing templates for GW detection, we will 
see how the ring-down phase could be matched to the inspiral phase in the EOB model, assuming that the QNMs are 
responsible for raising the frequency from the decoupling to the plateau. 

Finally, if we denote by merger the phase from roughly the decoupling point when M lo dec ^ 0.14-0.16 to the 
peak of Re[_2C22] when Mwpeak ~ 0.2 or the peak of radiation, the merger occurs in a very short time w 10-15M, 
corresponding to ~ 0.5-0.75 GW cycles. During this phase the frequency increases by « 45%, causing the GW 
spectrum to spread over a large frequency range (see Figs. and We shall discuss how this will affect the 

detectability of GWs from equal-mass binaries in Sec. IVIIII 

VII. EFFECTIVE-ONE-BODY APPROACH TO INSPIRAL-(PLUNGE)-MERGER-RING-DOWN 

The Taylor-expanded Hamiltonian for a two-body system was computed at 3PN order in Refs. [o^, [9l|. It took 
several years to compute the GW energy flux at 3.5PN order 25]. Before the 3.5PN dynamics was completed the 
Taylor-expanded FN predictions for the GW energy flux and the phasing of equal-mass binaries were not accurate 
enough to obtain robust predictions of the GW signal during the last stages of inspiral and plunge. For example, 
through 2.5PN order, the PN-approximants of some of the crucial ingredients entering the GW signal, such as the GW 
energy flux, differ significantly when evaluated at subsequent FN orders in the typical frequency band of ground-based 
detectors [33 ■ On the other hand, the signal-to-noise ratio of ground-based interferometers (especially the LIGOs) 
reaches its maximum around comparable-mass binaries. Thus, likely, the first detection may come from a coalescence 
of stellar-comparable-mass BHs merging in the most sensitive region of the detector's frequency band. 

At the beginning of 2000, in absence of NR results and under the urgency of providing templates to search for 
comparable-mass BHs, some resummation techniques of the FN series were proposed. In Ref. [3j|, the authors pro- 
posed the Fade resummation of the two-body energy and the GW energy flux, and in Refs. [l^l a specific resummation 
of the FN Hamiltonian was proposed, the so-called effective-one-body (EOB) Hamiltonian. Later on, in Ref. [28] the 
last stages of inspiral and plunge were modeled by combining the EOB Hamiltonian with the Fade resummation of 
radiation-reaction effects, providing the GW signal which includes effects beyond the adiabatic approximation. The 
EOB Hamiltonian was then extended at 3FN order without spin effects in Ref. [sol] and with spin effects in Ref. [3l| . 
More recently, the transition from inspiral to plunge including spin couplings has been modeled in Ref. [3^ . These 
analytical studies predicted that (i) the two-body motion would be quasi-circular throughout the last stages of inspiral 
and plunge, until the light-ring (see Fig. 1 in Ref. [1^), (ii) the ISCO for an equal-mass binary is a rather blurred 
concept, taking place roughly during half of a GW cycle (see Fig. 12 in Ref. [1^) and that (iii) the adiabatic plunge 
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FIG. 18: Features of the merger phase. We plot the frequency evaluated from the orbit and the wave, and the GW energy flux. 
We mark with circles the time when the common AH of the final BH first appears, when the binary separation reaches the 
light-ring of the final BH, the peak of the radiation flux (which occurs around 3 — 4Af before the peak in the amplitude of the 
waveform), when 50% of the energy and angular momentum have been radiated, and the time when 99% of the energy has been 
radiated (99% of the angular momentum appears to be radiated around 5 — IQM before this, though due to the oscillations in 
dJz/dt we are much less certain exactly when this occurs — see Fig |27p . The left panel refers to the d = 16 run and the right 
panel to the d = 19 run. 




FIG. 19: We compare the NR orbital frequency and the EOB orbital frequency obtained using the nominal x value from Table 
in The left panel refers to the run d — 16, the right panel to d = 19. The horizontal light continuous line marks the ISCO 
frequency, as predicted by the conservative 3PN-EOB Hamiltonian. The two dashed lines span the region in which a dynamical 
ISCO might be present. 
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FIG. 20: We compare the NR Re[-2C22] and the BOB Re[-2C22] waveforms obtained using the nominal x value from Table 
|I1 The left panel refers to the run d — 16, the right panel to d = 19. The two vertical dashed lines span the region in which a 
dynamical ISCO might be present. 



lasts only for almost one GW cycle (see Fig. 12 in Ref. [28jl. 

In Refs. [28l[3^. the authors also provided an example of the full waveform by modeling the merger as a very short 
(instantaneous) phase and by matching the natural end of the EOB plunge (around the light-ring) with the ring-down 
phase (see Ref. 92] where similar ideas subsequently developed also in NR). The matching was done using only the 
least damped QNM whose mass and spin were determined by the binary BH energy and angular momentum at the 
end of the EOB plunge. The choice of the light-ring at « 3Af for shifting the description between a (quasi-circular) 
binary motion and a deformed Kerr BH, was inspired by two considerations [28j|. First, in the test-mass limit, v <^ 1, 
Refs. [1^, [13] (see also Ref. [391]) realized a long time ago that the basic physical reason underlying the presence of a 
universal merger signal was that when a test particle falls below 3M (which is also the unstable light storage ring of 
Schwarzschild), the GW it generates is strongly filtered by the effecti ve p otential barrier centered around it. Secondly, 
for the equal-mass case v = 1/4, the close limit approximation [t^ [U [U, HI, HI] suggests a matching between the 
two-body and the perturbed-black-hole descriptions when the distance modulus /ip — 2, which would correspond to a 
Schwarzschild-like radial distance ~ 2.6M. 

For non-spinning, equal-mass binaries, the EOB approach predicts at 3PN order [3^: Ocnd/Afond — 0.77, Mend — 
0.974 M, an energy released of ~ 1.3%Af until t^nd and « 1% from the least-damped QNM phase [1^, Here 
and henceforth we denote by the subscript "end" the time at the end of the EOB plunge. Depending on the PN 
order the ending time can occur slightly before the conservative light-rin g. t c,„ci is the time at which the quasi-circular 
assumption used in building the EOB equations of motion breaks down [32'] . 

We show now some first-order comparisons between the numerical waveforms and the EOB waveforms, evaluated 
with the EOB conservative dynamics at 3PN order with spin-orbit and spin-spin effects included through 2PN order, 
and Pade radiation- reaction energy flux at 3.5PN order [12]. In contrast to the PN adiabatic model discussed in 
Sec. IIVB[ which should be used with confidence only until the last stable orbit, the EOB model extends beyond it 
through the plunge. 

In Sec. IIVBI a comparison of the PN orbital frequency with the numerical results was obtained by fitting to the 
orbital frequency lOc- In the EOB model the orbital frequency is obtained by solving the EOB equations of motion, 
thus it is more complicated to implement a least-square fit. Here, we simply determine by hand which initial EOB 
orbital frequency best matches, on average, the numerical orbital frequency ujc from the initial time to the light-ring 
or CAH (see Sec. IVip. and compute the corresponding wave. In Fig. [19] we show the results for the d = 16 and d = 19 
runs assuming the nominal x value from Table [H At the initial time t = 0, we find AI ujq = 0.047, Eq/M — 0.986, 
Jo/Af 2 ^ 0.896 for d = 16 and M cuq = 0.034, Eq/M = 0.988, Jo/M^ = 0.945 for d = 19. The evolution ends at 
Tend ^ 2.5 M, where Mwcnd ^ 0.16, E^nd/M ~ 0.971, Jcnd/M^ ~ 0.741 for the d ^ 16 run, and at r^nd ^ 2.5 M 
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FIG. 21: We compare the NR and EOB frequency and Re[_2C22] waveforms throughout the entire inspiral-merger-ring-down 
evolution. The data refers to the d — 16 run. 



where M Wond ^ 0.16, Eend/M ~ 0.971, Jend/M"^ ^ 0.737 for the d = 19 run. We notice that the EOB conservative 
ISCO for X = 0.063 is at Mwisco ~ 0. 096 [H i [Hi, rather close to the frequency range Mwdynisco = 0.078-0.097 
of the dynamical ISCO discussed in Sec. IIII D] In Fig. [201 we compare the NR and EOB Re[_2C22] waveforms. The 
two vertical dashed lines in Fig. [2D]mark the region during which a dynamical ISCO may be present. We compute 
the EOB waveform using Eq. ([^5)1 . which is valid in the adiabatic circular-orbit case. We checked that by relaxing 
this assumption and computing Re[_2C22] by taking derivatives of the binary quadrupole moment, the wave does not 
change much, except at the very end. 

By assuming the merger is a very short phase, the authors of Ref. [1^ simply joined the GW signal at the end of 
the inspiral to the least-damped QNM. As said above this modeling was inspired by the idea that once beyond the 
light-ring (i.e., inside the potential barrier), the GW emission is quickly dominated by the excitation of the QNM of 
the newly-formed BH. The choice of matching only one QNM inevitably creates a sudden jump of the GW frequency 
at the matching point. However, a smoother transition can be obtained by including higher overtones. As discussed 
in Sec. IVI[ the analysis done in Sec. IVl would suggest that the QNM production starts a bit later, around the peak of 
the radiation. At this stage we do not know whether a linear superposition of QNM can be responsible for raising the 
frequency from around the light ring (see Fig. llSp to the peak of radiation. In the spirit of an effective approach aimed 
at modeling the GW signal for detection, we push this idea further by including higher overtones when matching to the 
ring-down phase and discuss the consequences. The inclusion of higher overtones when matching to the ring-downn 
phase has also recently been adopted in Ref. fsBj, where the authors computed the transition inspiral-plunge-ring- 
down of a test particle in Schwarzschild. By including higher overtones, the authors could successfully match the 
exact numerical rise of frequency from the light ring to the least damped QNM as obtained from the Zerilli equation. 

First, we evaluate the BH mass and angular-momentum at the end of the EOB plunge, finding Mend = E^nd — 
0.971 M and acnd/Mend = Jcnd/E^^^ ~ 0.785. Then, we notice that those values are not the final BH mass and 
angular momentum because the binary has yet to emit energy and angular-momentum from the light-ring or CAH to 
the least-damped mode (see discussions around Fig. [TH]) . Future NR simulations will provide predictions for different 
mass ratios and spins. Here, guided by the results of Sec.lVIIand Fig.[T8l we assume that (Mend — M/)/Mond = 1.5% 
and (flcnd — a/)/acnd = 6%. Thus, we obtain Mj = 0.956 Af and af/Mj = 0.738. Using Ref. [76!|, we determine 
the frequency and the decay time of the fundamental mode and the first two overtones, finding Ma;220 = 0.576, 
Muj22i = 0.565, Mw222 = 0.545, T220/M = 0.0828, r22i/A/ = 0.250 and T222/M = 0.422. We then determine the 
three unknown amplitudes and three unknown phases of the three QNMs by imposing the continuity of the frequency 
^D2, Eq- (|23p . the wave and its first five derivatives at t = tend- We note that this matching procedure is rather 
sensitive to the time of matching because the frequency is increasing very quickly around tend- In Fig. [^l] we compare 
the frequency and the inspiral-(plunge)-merger-ring-down wave Re[_2C22] of the EOB model with the NR results for 
the case d = 16. The EOB ring-down frequency is computed from Eq. where we used in Im[_2C22] the same three 
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amplitudes and three phases of Re[_2C22]- Were we to use in Ini[_2C22] the three ampHtudes and phases obtained 
by matching it at tondj we would not obtain a good result. This is due to numerical errors introduced by matching 
separately Re[_2C'22] and Im[_2C'22] at iond- 

As seen in Fig. [511 by matching the fundamental QNM and the first two overtones, the frequency transition becomes 
smoother, but nevertheless it differs from the NR frequency tox. As we shall see in the next section, this effective way 
of including a short in time, but spread in frequency, merger phase, can mimic the frequency spread of the power 
spectrum of the NR waves, though with a slightly different power law. 

VIII. DETECTABILITY OF THE SIGNAL 

In this section we compute the Fourier transform of the numerical waveforms, compare the results with the analytical 
predictions and give an estimate of the optimal signal-to-noise ratio (SNR) for ground-based and space-based detectors. 

Frequency domain PN templates for the inspiral phase are generally computed in the so-called Stationary Phase 
Approximation (SPA). They read 

'.spaW^ .4 ^ = ^^ (45) 

where / is the frequency of the GWs, A4c = ly^^^Ad is the chirp mass and Z?l is the luminosity distance to the source. 
In Eq. (|45p we have adopted the standard "restricted PN approximation" , in which the amplitude is expressed to the 
leading order in a PN expansion while the phasing ^spa(/), is expressed to the highest PN order available. The phase 
is currently known through 3.5PN order. Here, we are only interested in computing the amplitude of the Fourier 
transform of the signal and investigating how and when it starts deviating from the Newtonian prediction /~^/®. We 
shall analyze the comparisons between the Fourier transform phases in the future. 

In the left panel of Fig. [21] we plot the Fourier transform amplitudes of the numerical waveform for the three runs, 
for a (15 + 15)Mq binary which is a typical source for LIGO/VIRGO/GEO/TAMA. Using Table I we find that for 
this binary mass the initial GW frequency in the three runs is 121 Hz, 89 Hz and 71 Hz (vertical dot lines in the 
left panel of Fig. I22p . To compute the Fourier transform we extrapolate the numerical waveforms at earlier time, for 
almost 6 X 10^ m, by attaching to it the 3PN-adiabatic model which best-fits it. We compute the Fourier transform in 
three different ways, from Re[^'4] and /i+ extracted along the direction perpendicular to the orbital plane, and from 
-2(^22- Besides a normalization factor, the amplitudes computed from Re[^'4] and -2^22 must agree in the inspiral 
phase, where they satisfy more and more the restricted PN approximation, but they can differ in the last part of 
inspiral, merger and ring-down, when non- linear effects and higher harmonics can become important. 

From Fig. I22[ we see that at low-frequency, during the last stages of inspiral, the amplitude can be approximated 
by the Newtonian amplitude Z^^^^, with small bumps maybe due to the presence of eccentricity. At higher frequency, 
during the merger and ring-down the slope changes to /~". At this stage we cannot uniquely determine the slope 
index or say if there is more than one change in the slope. We estimate n « 0.6-0.8 ^ and notice that the change in the 
slope can occur as early as the beginning of the plunge AIuj ~ 0.1 or as late as the decoupling time M uj ^ 0.16. As an 
example in Fig. [Hjwe show the case n — 2/3. Finally, the signal drops at higher frequencies, around the frequency of 
the fundamental QNM, i.e., 620 Hz for the (15 -I- 15)Af0. Even if the merger occurs in a short time, corresponding to 
~ 0.5-0.75 GW cycle, the frequency increases very quickly during this phase and then approaches the frequency of the 
fundamental QNM. As a consequence, the Fourier transform signal spreads over a large frequency band 300-600 Hz. 
We also computed the Fourier transform amplitude from ^'4 extracted along directions 6^0, and found differences 
in the merger-ring-down amplitude slope. We shall discuss those interesting features in a future publication. 

In the right panel of Fig. [HI we show the sky averaged SNRs versus total mass, for an equal mass binary at 100 Mpc, 
and for initial LIGO. The dashed-dot and dashed curves are computed assuming the SPA inspiral signal (jlS)) until 
the Schwarzschild ISCO /isco = 4400/ (M/Mq) and the ICO predicted by adiabatic PN theory at 3PN order 
respectively. The average SNR for one detector for an SPA signal is computed using 

The continuous light curve in Fig. [22| is computed from the EOB inspiral-(plunge)-merger-ring-down wave shown 
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Similar results were also obtained independently in Ref. (86|] 
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FIG. 22: In the left panel we plot the amplitude of the Fourier transform of the numerical waveform for a binary with (redshift) 
mass (15 + 15)Mq. The lower horizontal axis marks the gravitational wave frequency in Hz, while the upper axis marks the 
dimensionless orbital angular velocity that coincides with the instantaneous GW frequency. The vertical dotted lines mark the 
frequencies at which the runs start. In the right panel, we show the average SNR for one detector versus the (redshift) total 
mass for an equal-mass binary at 100 Mpc. 



in Fig. The change in the slope that we observe in the left panel of Fi g. toget her with the inclusion of the 
signal beyond the ISCO, causes an increase of the SNR for large masses (28l l3ll |3^ |93|. For total masses lower than 
(15 + 15)Mq, the end of the inspiral occurs around the most sensitive LIGO frequency, while the merger-ring-down 
is pushed to higher frequency where the sensitivity is much lower. As a consequence, the average SNR which includes 
merger and ring-down phases does not differ much from the one that includes only the inspiral phase. Astrophysical 
observations and theoretical predictions suggest that stellar mass BHs have a total mass ranging between 6-30 Af© . If 
binary BHs of larger total mass exist, they could be detected by initial LIGO with very high SNR. The FOB model 
with the instantaneous matching to three QNMs predicts a SNR close to the NR result. It is smaller because the 
effective matching of the inspiral to the three QNMs (see Fig. [^Tjl gives a Fourier transform amplitude which extends 
until the QNM frequency, but with a slope index slightly larger than « —0.6-0.8. 

In the left panel of Fig.[23]we plot the amplitude of the Fourier transform for a signal typical of LISA, a (1O^-|-1O^)M0 
supermassive BH binary. In the right panel we show the average SNR for one Michelson LISA configuration versus 
the total (redshift) mass for an equal-mass binary at 3 Gpc. The dip in the plot is due to the WD- WD confusion 
noise [87]. Due to the inclusion of merger and ring-down phases, the SNR increases considerably for total masses 
larger than 2 x IO^Mq. We notice that Fig. [23] is consistent with Fig. 7 of Ref. [76] where the authors computed the 
SNR due to the ring-down phase, assuming ^ 3%to of energy released during the merger. 

IX. CONCLUSIONS 

In this paper we have analyzed the data of several numerical simulations of the inspiral-merger-ring-down of an 
equal-mass binary carrying a very small spin aligned with the orbital angular momentum. 

The combination of several effects including (i) limited resolution, (ii) relatively close initial configurations and (iii) 
lack of diagnostics to measure and compensate for possible coordinate artifacts, make it impossible to claim very 
high accuracy in the comparisons with analytical models and the analysis of the merger waveform. Nevertheless, the 
resolution studies performed in Sec. IIIII suggest that we are in the convergent regime, and so meaningful conclusions 
can be drawn from the data. Furthermore, the consistency with which several quantities have been measured by 
independent means suggest that adverse gauge effects are minor and will not affect many of the conclusions reached. 
In particular, as Tables [TTTl and fVl show . the final mass and angular momentum extracted from the ring-down are very 
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FIG. 23: In the left panel we plot the amplitude of the Fourier transform of the numerical waveform for a binary with (redshift) 
mass (10^ + 1O®)M0. The lower horizontal axis marks the gravitational wave frequency in Hz, while the upper axis marks the 
dimensionless orbital angular velocity that coincides with the instantaneous GW frequency. The vertical dotted line marks the 
frequency at which the d = 19 run starts. In the right panel we show the average SNR versus the (redshift) total mass for an 
equal-mass binary at 3 Gpc {z = 0.54). 



close to their values measured through AH properties. Also, the orbital frequency measured via AH motion is close 
to the one extracted from the GW, as Fig. [7] shows. 

We found that one of the dominant forms of numerical error is a slow drift in the phase of the waveform. With 
respect to detectability in a GW burst-search this error does not seem to be very significant, as the cumulated phase 
error up to the peak in the radiation can be factored out by a constant phase shift and the subsequent coalescence/ring- 
down waveform does not seem to be significantly affected by prior phase error. However, for matching to FN models 
at earlier time and parameter estimation this error is significant, and directly translates into uncertainties in matching 
parameters. 

As Figs. [5] and [3] show, the numerical evolution is characterized by a strong initial pulse of radiation. The timing 
of the pulse suggests that it is associated with the assumption of conformal flatness. In fact it has very similar 
characteristics to the initial pulse seen in scalar field collapse generated binaries [1, [sHi, which also begin with a 
conformally flat spatial metric. Fortunately, the effects of this pulse of radiation seem to diminish rapidly. The initial 
data does not put the binary on a dean circular inspiral path. As discussed in Sec. IIIIDI (see Fig. [51), the trajectory 
clearly oscillates about the desired trajectory. The effects of this oscillation can be interpreted as a small eccentricity 
of the initial orbit. We estimated it to be e '--^ 0.02 for the d = 19 case. While these oscillations could be due to an 
actual eccentricity in the initial data, they could also be due to a lack of an appropriate initial radial momentum. The 
current simulations cannot determine which effect, if either, is most significant. However, while there are oscillations 
in the inspiral trajectory, the resulting dynamics is still adiabatic (e.g., see Fig. [5]) and modeled well by circular orbits. 

Concerning the inspiral phase, as described in Sec. IIV Al it is remarkable how well the Newtonian quasi-circular 
approximation can match the numerical signal of an equal-mass binary. For this we mean that if the orbital phase 
is modeled well, the leading order Newtonian term in the expansion of the waveform is able to model both the 
amplitude and phase of the GW quite accurately until close to the time of merger. It is also striking how, despite 
possible coordinate artifacts, a GW signal computed with the leading quadrupole formula using the coordinate motion 
of the AH's matches the numerically extracted wave to a reasonable degree (see Fig. [3]). 

As discussed in Sec. IIVBI the FN adiabatic phase and frequency matches ujc (frequency computed from the AH 
motion) well. In particular, if the analytical expressions for oj ((28)) and (j) (|3T|) are used, the 3PN-approximant best fits 
the data, whereas if the expression for Co p2p is solved numerically, a 3.5PN-approximant also matchs the numerical 
data to a similar level. The accuracy of the PN-adiabatic model improves with increasing binary separation. We 
expect it deviates from non-adiabatic models when approaching the last stable orbit. Nevertheless, we found that if 
the PN-adiabatic model is extrapolated even to the formation of the CAH, it gives reasonable results. This is due to 
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the fact that the numerical plu nge cycle is very short and still quasi-circular for an equal-mass binary, as originally 
predicted in the EOB model [2^. In Sec. lVIII we compared the numerical results with the EOB model. Again we found 
that the 3.5PN-approximant best matches the data. The blurred IS CO p hase of half of a GW cycle and subsequent 
plunge phase of almost one GW cycle predicted in the EOB model [28] seems to be consistent with the numerical 
results. By properly adjusting the mass and angular-momentum of the BH at the end of the EOB plunge, to account 
for the energy and angular-momentum released during the merger, and by including three QNMs, we could match 
the EOB inspiral- (plunge) wave to the ring-down wave (see Fig. I2ip . This matching procedure is very sensitive to 
the time of matching and it is only partially effective (see the difference in the GW frequency in Fig. [^T|l since it does 
not capture the details of the merger, but could be improved in the future. 

As Fig. [7] shows, decouples from the dominant GW frequency just before the point of formation of the common 
AH, around the light ring. We conjecture that this decoupling time marks the transition point between inspiral- 
(plunge) and merger. We found that the merger accounts for only a brief time 10-15M) compared to the inspiral 
(arbitrarily long) or ring down phases (~ 30M as measured by adding T220 to the time at which higher-order modes 
and overtones become insignificant). However, the dominant GW frequency rises very quickly and spans a significant 
range of frequencies during the merger, as can be seen in the left panel of Fig. [22l With current data, the peak in 
1-2(^22! (which occurs a few M after the peak of the radiation) seems to be a natural point to mark the transition 
between the merger and ring-down phases. It is possible that the higher-order ring-down modes/overtones (with 
lower frequencies than 0^220) are excited by resonance with the dominant GW frequency as it rises during the merger. 
An open and important question remains whether it would be possible to model the merger by using some kind of 
non-linear modification of the onset of each QNM. 

Concerning the ring-down phase, we extracted the fundamental QNM and first few overtones. Results are shown 
in several Tables of Appendix ID] and Figs. [T51 and [TBI We fit to the individual _2Cfm modes instead of to ^4 directly 
because fitting to _2C^m includes information from all directions. The fit to each _2Cfm can, in principle, include 
QNMs for all values of ^ > |to| but only the negative frequency modes with azimuthal index m and positive frequency 
modes with index —m. For -2^22, we found only the £ = 2,m — 2 QNMs are significant. As seen in Fig. \17\ we 
could fit the ring-down signal to times slightly before the peak in |_2C22| by using n = . . . 3. For _2C32, we found 
that we must use both the £ — 2,m — 2 and £ — 3,m — 2 QNMs, but other modes do not contribute at a significant 
level. Also in this case we can again fit the ring-down signal to times slightly before the peak in |_2C32| by using 
n = ... 2. More sophisticated ring-down fitting techniques will be helpful in gaining a better understanding of the 
transition from merger to ring down. 

Quite interestingly throughout the inspiral-merger-ring-down the balance equation dE/dt = LodJ/dt holds quite 
well on average, as Fig. [27] in Appendix [BI shows . with a single frequency always dominating the entire evolution (see 
Fig. [ZD . 

The analysis of the inspiral-merger-ring-down suggests that it should be possible to come up with good hybrid 
numerical/analytical waveforms, or even complete analytical waveforms where the full numerics guides how we need 
to patch the inspiral and ring-down waveforms together. To this end, it will be very important to devise a simple 
model of how the QNM's are excited during the transition regime. Of course, all of this will be moot if the relative 
simplicity of this merger scenario breaks down for more interesting initial conditions (mi 7^ m2, 6*1 7^ 5'2 and 
non-aligned with the orbital angular- momentum) , so it will be important to numerically evolve more varied classes of 
initial data. 

Finally, during the relatively brief merger phase the dominant GW frequency rises quickly, generating a signal in the 
Fourier domain that is rather spread out in frequency. The left panels in Figs. [22| and [2S| indicate a change of slope in 
the signal Fourier amplitude h. The slope during inspiral is —7/6. The slope during the transition merger-ring-down 
seems to be ~ —0.6-0.8. The inclusion of the merger-ring-down signal increases the SNR for large binary masses. If 
binary BHs of mass larger than 40Mq exist, they could be detected by one single LIGO with SNR up to 15 at 100 
Mpc. LISA could detect supermassive BHs of masses 2 x 10^-10^, with SNR up to IC* at IGpc. 
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FIG. 24: Magnitudes of several sub-dominate components of the d = 19 waveform relative to the dominant 1 — 2, |m| = 2 
component. Modes up to ^ = 6, \m\ = 6 were extracted — those not present in the plot were several times (at least) smaller 
than any of the modes shown. 



APPENDIX A: ADDITIONAL MULTIPOLE MOMENTS 



For the majority of comparisons in this paper we focused on the dominant £ = \m\ = 2 spherical harmonic 
components of a waveform, as for quasi-circular, equal mass binaries that are initially co-rotating one does not expect 
other multipoles to be present to any significant degree. In Fig. [23]we plot the magnitudes of several higher multiple 
moments of the waveform relative to the -2C2.2 multipole moment for the d — 19 h/2 run (the d — 13 and d — 16 have 
a similar spectrum). The next-to-leading order component of the waveform is the _2C4^4 moment, which is about 
an order of magnitude smaller than -2^2, 2- At up 5 times smaller than -2(^4,4 are the -2Ce>3,2, -2(^6,4 and _2C6,6 
modes. Interestingly, during the ring-down part of the waveform the relative magnitude of the -2C'f>3.2 modes grow 
and become almost as significant as the -204^4 mode. We only calculated components of \E'4 up to ^ = |m| — 6 — all 
other modes not shown (including the axisymmetric m = modes) are at least of factor of 3 times smaller than any 
of the modes in the plot. Note however that the sub-dominant components of the waveform are more susceptible to 
gauge effects and numerical truncation error in the solution as their relative amplitudes are so small; hence one should 
be wary of drawing any significant conclusions from Fig. [24l 

Even though the higher multipole moments are small, it will still be interesting to compare them to the predictions 
of perturbative calculations. Furthermore, for unequal mass binaries and/or binaries with significant initial spins — 
in particular with spin components that arc not aligned with the orbital angular momentum — the higher multiple 
moments will play a much more important role in the description of the waveforms. We leave it to future studies to 
analyze these modes in more detail. 



APPENDIX B: ENERGY AND ANGULAR MOMENTUM FLUX 



The GW energy flux emitted by a binary moving along an adiabatic sequence of circular orbits is currently known 
through 3.5PN order [2^ for non-spinning BHs and through 2.5PN for spinning BHs 70]. It reads 

= |,2(Mc.)W3|l+(^_^_|,^(M^)2/3 + 4^(^^^) 

44711 9271 65 oA/,. x4/3 f 8191 583 \ 
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FIG. 25: In the left plot we compare the GW energy fluxes of the numerical simulation and PN-adiabatic model when the 
best-fit for the orbital frequency is used (see Sec. lIVB|l . In the right plot, we show the same comparisons but we divide the fluxes 
by the Newtonian GW flux. For the numerical case we show the results assuming circular orbits, i.e., -Fkcwt = 32/5 v^^Mlli)^^^^ 
and not assuming circular orbits, i.e., -Fisiewt ~ S2/5u^{Mui)^ (M/r)~*. The data refers to the d = 16 run. 
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The numerical energy flux is computed in terms of the mode amplitudes _2Cfm(i) as 
where Dimit) is a dimensionless first time integral of -^2Cfm(t) defined by 



(Bl) 



(B2) 



M 



dt' -2Cim{t')- 



(B3) 



In the left panels of Figs.[25l[26l we compare the PN-adiabatic and numerical fluxes for the runs d = 16 and d = 19. 
We compute the PN-adiabatic flux from Eq. (|Bip . using the best-fit w derived in Sec. IIV Bl To better pinpoint the 
differences, we show in the right panels of Figs. [25l [26l the GW energy fluxes normalized to the Newtonian GW flux, 
assuming the Keplerian relation uP'r' jM = 1 valid for circular orbits. In addition, we show a curve where the flux has 
been calculated without assuming the Keplerian relation. We notice that for the case d = 19, except for the very first 
part of the evolution, the PN-adiabatic flux averages the numerical one until ~ 30-50M before the CAH forms. For 
the case d = 16, the PN-adiabatic flux computed from the best-flt PN-adiabatic u! always overestimates the numerical 
one. 
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FIG. 26: In the left plot we compare the GW energy fluxes of the numerical simulation and PN-adiabatic model when the 
best-fit for the orbital frequency is used (see Sec.|lVB}. In the right plot, we show the same comparisons but we divide the fluxes 
by the Newtonian GW flux. For the numerical case we show the results assuming circular orbits, i.e., Fncwi = 32/5 ^^{Mu 



and not assuming circular orbits, i.e., -Fncwi ~ 32/5 u (Muj) (M/r) 
dynamical ISCO could exist. The data refers to the d = 19 run. 
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Figures 1261 show that the frequency obtained by fitting the numerical frequency, as done in Sec. lIVBl not only 
provides a GW signal that matches the numerical one quite well, but also provides an energy flux that is consistent 
with the numerical energy fiux. However, this study does not give hints on which analytical model and/or PN order 
best fits the numerical flux. We leave for future work a more detail study of the comparison between numerical and 
analytical energy and angular-momentum fluxes, which includes fully relaxing the assumption of circular orbits and 
comparisons with the Pade resummed fluxes. 

In PN theory it is possible to show that the following relationship holds between the radiated angular momentum 
flux dJz/dt and energy flux dE/dt for circular orbits: 

dE dJz . X 

— = w — ^ . B4) 

dt dt ^ ' 

The above equation is used in building the EOB equations of motion [28l . Is^ l . The numerical angular momentum flux 
is computed in terms of the mode amplitudes -2Cim{t) as 

^ - ^^J2'^(.MDUt)EUt)]), (B5) 

im 

where Eim{t) is a dimensionless second time integral of -~2Cim{t) defined by 

1 /■* 

El„,{t) = — / dtDirnit'). (B6) 



MJo 

In Figure [27| we plot these two quantities calculated from the angular momentum and energy fluxes of the numerical 
simulation, using the frequency lu]J2 P3|) . The oscillations we see in the numerical dJz/dt are, we believe, due in 
part to improper initial conditions in the double-time integral of 'I'4 used to obtain dJz / dt (we used initial conditions 
that are valid at t = — oo); however, on average (|B4I) appears to hold to a reasonable degree throughout the entire 
evolution. 
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APPENDIX C: EFFECT OF EXTRACTION RADIUS ON MEASURED WAVEFORM 

There are many coordinate related issues in extracting GWs from the simulation, most of which have not been 
rigorously addressed for this set of runs. The issues can be summarized by the question how does the gauge used in 
the simulation affect the assumed relationship between ^'4 and the gravitational wave strain given in Eq. |_?^| ) ?. Most 
certainly there are differences that decay by some power of r, as only in the limit as r — s- 00 is (jl4p strictly satisfied. 
To test whether such differences are present and estimate how significant they are, we can examine rM"^4^{t, r, 6, 0), 
or equivalently _2Cfm(i, r) as defined in Eq. (fTT]) . as a function of extraction radius r. In the wave zone of a "good" 
coordinate system -2Cim{t,r) = -2Ctm{t — r) for large r, and so by comparing the waveform at several extraction 
radii will give some hints as to the adequacy of the coordinates. 

Figures [28H301 show two plots each of -2^2.2(^,7") at four extraction radii, r — I2.5Af, 25M, 37.5Af and 50M, for 
each of the three simulations (h/2 resolution). The plots on the left show the early part of the waveform, shifted by 
an appropriate amount in time assuming that the wave is propagating with unit velocity. As can be seen, going from 
an extraction radius of 37. 5M to 50M the shifted waveforms all overlap quite closely, in consistency with wave-like 
propagation. However, these time-shifts do not give such a good match in the part of the waveform associated with 
the coalescence and ring-down of the binary. This is evident from the plots on the right, where now the time shift has 
been chosen to give the best possible phase overlap around peak amplitude. What these latter plots suggest is that 
the gauge changes by a small amount with time in the wave extraction regime of the simulation, and the amount of 
change is apparently correlated with the amplitude of the GWs being emitted. 

These effects are summarized in Table IVH where we give the average coordinate propagation speeds between the 
different extraction radii during the inspiral and coalescence/merger portion of the signal. Another apparent gauge 
pathology seen in Figs. [28l - [30l is that near the peak amplitude the assumed 1/r decay in the wave amplitude [factored 
out in Eq. (|f 7p ] does not describe the situation as well as in the inspiral regime. This could happen (for example) 
if the coordinate sphere r = begins to deviate from a geometric sphere of radius r^; such a change in gauge could 
also affect the coordinate velocity of the wave, and in fact the trend of decreasing amplitude with extraction radius 
is consistent with the decrease in wave speed (though cannot by itself account for all the change in velocity). 



APPENDIX D: TABLES 



This appendix contains tables from the QNM ring-down analysis presented in Sec.|Vl Tables FVIIl I VIII and IIXI are fit 
parameters for the -2C22 component of the waveform from the d = 13, 16 and 19 cases respectively; similarly Tables 
1X1 1x1 and IXIII are fit parameters for the _2C32 component of the waveform for the same set of separations. 
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FIG. 28: A component of the waveform of the d=13 (h/2) run measured at different extraction radii, and shifted in time to 
account for wave propagation speed. The figure to the left has time shifts assuming unit coordinate speed, which for large 
extraction radii seem to be a good assumption for the inspiral part of the wave though produce some mismatch around peak 
amplitude; the figure to the right depicts the time shift necessary to produce the best overlap about peak amplitude. See Fig's 
1291 and [30l for similar plots for d=16 and d=19 respectively. 




FIG. 29: Time-shifted waveforms versus extraction radius for d=16 (see the caption describing the the corresponding d — 13 
plot in Figl2HI 
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TABLE VI: The average GW propagation speed between various extraction radii for the inspiral part of the wave and the 
coalescence/ring-down ("peak") part — see Fie:.'s [28l tol30l 
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FIG. 30: Time-shifted waveforms versus extraction radius for d=19 (see the caption describing the the corresponding d = 13 
plot in Fig l28|l 



APPENDIX E: RESULTS FOR d = 13 

This appendix contains some figures showing results from the d = 13 run that were excluded from the main text for 
brevity. Figure [311 shows orbital angular frequencies extracted using several methods (Sec. IIV A| . and a comparison 
of the numerical and NQC inspiral waveforms (Sec. IIV A|) . Fig. [32] shows a comparison of frequencies and waveform 
components between NR and various analytical counterparts (Sec. IIVB|) . Fig. [551 shows results from QNM extraction 
(Sec. El, and Fig. [34] plots the dominant frequencies during the ring-down phase (Sec. EJ and identifies features of 
the merger (Sec. IVI|) . 



[1] A. Abramovici et al., Science 256, 325 (1992); [http: //www. ligo . caltech. eduj 
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TABLE VII: Fit parameters for -2C22 for the d = 13 case, trim denotes the time a/ter the peak in I-2C22I at which the fitting 
started. 
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N = 2 
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TABLE VIII: Fit parameters for -2C22 for the d = 16 case, tr/m denotes the time after the peak in I-2C22I at which the fitting 
started. 
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TABLE IX: Fit parameters for -2C22 for the d = 19 case. tr/M denotes the time after the peak in |_2C22| at which the fitting 
started. 
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TABLE X: Fit parameters for -2C32 for the d = 13 case. tr/M denotes the time after the peak in I-2C22I at which the fitting 
started. 
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TABLE XI; Fit parameters for -2C32 for the d = 16 case, tr/m denotes the time after the peak in I-2C22I at which the fitting 
started. 
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TABLE XII; Fit parameters for -2C32 for the d = 19 case. tr/M denotes the time after the peak in j-2C22l at which the fitting 
started. 
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FIG. 31: (left panel) Orbital angular frequencies for the d = 13 case. See Fig.[7]for a full description, (right panel) Comparison 
of numerical and NQC inspiral waveforms for the d = 13. See Fig. [8] for a full description. 
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FIG. 32: (left panel) Comparison of the NR and three analytical orbital frequencies — see Fig. [10] for a full description, (right 
panel) Comparison of the NR and analytical Re[--2C22] — see Fig. 1111 for a full description. 
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FIG. 33: The left panel shows a comparison of numerical and QNM -2C22 ring-down waveforms for d = 13 — see Fig. 1151 for 
a full description. The right panel shows a comparison of numerical and QNM -2C32 ring-down waveforms for d = 13 — see 
Fig. US] for a full description. 
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FIG. 34: (left panel) Dominant frequencies during the ring down for the d = 13 case evaluated using several methods. See 
Fig. [13 for a full description, (right panel) Features of the d = 13 merger phase. See Fig. [TS] for a full description. 



